Key Results to Accompany Approach-Avoid RLWM Paper
Peter Frank Hitchcock - Translational Lab, Emory University
2024-10-11
- Set up packages, function sourcing, and plotting parameters
- Behavioral data
- Figure 2 and associated stats — correct (ie. reward) learning curves and capture by model, RT, and train-test correlation by set size
- Test correct (ie. reward)
- Incorrect trials
- Individual differences - depression/anxiety and rumination
- Comparison to choice kernel rather than RL model
Set up packages, function sourcing, and plotting parameters
sapply(c(
"rjson",
"data.table",
"dplyr",
"ggplot2",
"stringr",
"purrr",
"foreach",
"patchwork",
"testit",
"lme4",
"lmerTest",
"latex2exp",
"brms"#,
#"tidybayes" # Patchwork/tidybayes ggplot dependency conflict, so source tidybayes directly
),
require, character=TRUE)## Loading required package: rjson## Loading required package: data.table## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: ggplot2## Loading required package: stringr## Loading required package: purrr##
## Attaching package: 'purrr'## The following object is masked from 'package:data.table':
##
## transpose## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: patchwork## Loading required package: testit## Loading required package: lme4## Loading required package: Matrix## Loading required package: lmerTest##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step## Loading required package: latex2exp## Loading required package: brms## Loading required package: Rcpp## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:lme4':
##
## ngrps## The following object is masked from 'package:stats':
##
## ar## rjson data.table dplyr ggplot2 stringr purrr foreach
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## patchwork testit lme4 lmerTest latex2exp brms
## TRUE TRUE TRUE TRUE TRUE TRUEsf <- function() invisible(sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source)) # Source all fxs
sf()
DefPlotPars()Behavioral data
packageVersion("brms")## [1] '2.21.0'Note: These and other code paths need to be adjusted / dir
structure needs to be reconstructed for public data. The public data
also excludes demographics and, due to file size, simulations
(simulations can be rerun using s.R)
Read in learn and test df
learn_df <- read.csv("../data/learn_df.csv")
test_df <- read.csv("../data/test_df.csv")demogs <- read.csv("../data/demogs_deident.csv")Demographics
summary(demogs$Age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 19.00 30.00 37.00 38.42 47.00 65.00mean(demogs$Age); sd(demogs$Age)## [1] 38.42029## [1] 11.42298table(demogs$Sex)##
## DATA_EXPIRED Female Male Prefer not to say
## 1 139 133 3cat("\n Percent female, male, expired, Prefer not to say: \n",
139/nrow(demogs),
133/nrow(demogs),
1/nrow(demogs),
3/nrow(demogs))##
## Percent female, male, expired, Prefer not to say:
## 0.5036232 0.4818841 0.003623188 0.01086957table(demogs$Ethnicity.simplified)##
## Asian Black DATA_EXPIRED Mixed Other White
## 12 26 4 19 12 203cat("\n % Asian, Black, expired, Mixed, Other, White: \n",
12/nrow(demogs),
26/nrow(demogs),
4/nrow(demogs),
19/nrow(demogs),
12/nrow(demogs),
203/nrow(demogs)
)##
## % Asian, Black, expired, Mixed, Other, White:
## 0.04347826 0.0942029 0.01449275 0.06884058 0.04347826 0.7355072Mean delay by set size
learn_df %>% group_by(set_size) %>% summarize(mean(delay, na.rm=TRUE))## # A tibble: 5 × 2
## set_size `mean(delay, na.rm = TRUE)`
## <int> <dbl>
## 1 1 0
## 2 2 0.832
## 3 3 1.67
## 4 4 2.48
## 5 5 3.33Load parameter results from best model
## Write into one file and confirm 8/1/24 bug fix was minor
# 11325 and 15031 are the 8/1/24 runs with the minor bug/illogical scaling of epsilon correct
# m35_v1 <- read.csv("../model_res/opt/best/BEST__RunRLWMPRew11325.csv")
# m35_v2 <- read.csv("../model_res/opt/best/BEST__RunRLWMPRew15031.csv")
#ComparePars(m35_v1$nll, m35_v2$nll) # confirm essentially identical
#m35 <- rbind(m35_v1, m35_v2) %>% group_by(ID) %>% slice(which.min(nll))
# And basically identical to older version confirming was a minor bug
# m35_old <- read.csv("../model_res/opt/best/before_8-1-reruns/BEST_m35_RunRLWMPRew.csv")
# ComparePars(m35$nll, m35_old$nll) # confirm essentially identical
#write.csv(m35, "../model_res/opt/best/BEST__m35_RunRLWMPRew-8-1-24-epsilonfixed.csv")m35 <- read.csv("../model_res/opt/best/BEST__m35_RunRLWMPRew-8-1-24-epsilonfixed.csv") # the final model with the bug/nonsensical implementation of epsilon fixed on 8/1learn_m <- learn_df %>% group_by(ID) %>% summarize(m=mean(correct))
test_m <- test_df %>% group_by(ID) %>% summarize(m=mean(correct))
assert(all(learn_m$ID==m35$ID))
assert(all(test_m$ID==m35$ID))And simulations from model
m35_s <-
read.csv("../model_res/sims/SIM_RunRLWMPRew53679.csv")N simulations
length(unique(m35_s$iter)) ## [1] 50m35_s_learn <- m35_s %>% filter(type=="learning")
m35_s_test <- m35_s %>% filter(type=="test")Figure 2 and associated stats — correct (ie. reward) learning curves and capture by model, RT, and train-test correlation by set size
Plot the empirical proportion correct and its capture by the model
pcor_ss <- data.frame(learn_df %>%
group_by(set_size, stim_iter, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size', 'stim_iter'. You can override
## using the `.groups` argument.pcor_ss_err <- Rmisc::summarySEwithin(pcor_ss,
measurevar = "m",
withinvars = c("set_size", "stim_iter"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_size, stim_iterpcor_ss_m <- data.frame(learn_df %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.emp_p1 <- ggplot(pcor_ss_err, aes(x=stim_iter, y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") + tol +
ggtitle("Empirical") + tp## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Confirm summarySEwithin m and un-normed one are essentially the same
pcor_ss_m <- data.frame(learn_df %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.assert(all(round(pcor_ss_err$m, 4)==round(pcor_ss_m$m, 4)))pcor_ss_sim_m35 <- data.frame(m35_s_learn %>% group_by(set_size, stim_iter) %>%
summarize(m=mean(corrects), n()))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_ss_sim_m35_iter <- data.frame(m35_s_learn %>% filter(iter %in% c(1:30)) %>% group_by(stim_iter, set_size, iter) %>%
summarize(m=mean(corrects), n()))## `summarise()` has grouped output by 'stim_iter', 'set_size'. You can override
## using the `.groups` argument.sim_m35_p1 <-
ggplot(pcor_ss_sim_m35, aes(x=as.factor(stim_iter), y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
geom_jitter(data=pcor_ss_sim_m35_iter, aes(fill=as.factor(set_size)), color="black", height=0, width=.2, alpha=1, size=2, pch=21) +
geom_point(aes(fill=as.factor(set_size)), color="black", size=6, pch=21, alpha=.7) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("") +
tp + ggtitle("Simulated")emp_sim_perf <- emp_p1 + sim_m35_p1 emp_sim_perf
#ggsave("../paper/figs/pieces/fig2_emp_sim_perf.png", emp_sim_perf, height = 3.5, width=11, dpi=300)Before/after break regressor
end_p1_begin_p2_df <- learn_df %>% filter(stim_iter %in% c(5, 6))
#end_p1_begin_p2_df$before_after_break <- end_p1_begin_p2_df$stim_iter
end_p1_begin_p2_df[end_p1_begin_p2_df$stim_iter==5, "before_after_break"] <- -1
end_p1_begin_p2_df[end_p1_begin_p2_df$stim_iter==6, "before_after_break"] <- 1Statistical models of performance effect
Main effects of set size and stim iter and interaction reflecting parametric effect over time in both phase 1 and phase 2
summary(pcor_phase_1_v2 <-
glmer(correct ~ scale(set_size)*scale(stim_iter) +
(scale(set_size)*scale(stim_iter)|ID),
data=learn_df %>% filter(phase==1),
family="binomial", control = glmerControl(optimizer = "bobyqa"))) ## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(set_size) * scale(stim_iter) + (scale(set_size) *
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 1)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 44084.9 44205.7 -22028.5 44056.9 41236
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.8169 -0.7807 0.3279 0.6811 1.7107
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.19502 0.4416
## scale(set_size) 0.04951 0.2225 -0.08
## scale(stim_iter) 0.10059 0.3172 0.90 -0.23
## scale(set_size):scale(stim_iter) 0.02250 0.1500 -0.16 0.69 -0.53
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.99731 0.03038 32.83 <2e-16 ***
## scale(set_size) -0.54013 0.02075 -26.03 <2e-16 ***
## scale(stim_iter) 1.07075 0.02424 44.18 <2e-16 ***
## scale(set_size):scale(stim_iter) -0.37241 0.01824 -20.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) -0.201
## scl(stm_tr) 0.774 -0.297
## scl(s_):(_) -0.228 0.587 -0.417summary(pcor_phase_2_v2 <-
glmer(correct ~ scale(set_size)*scale(stim_iter) +
(scale(set_size)*scale(stim_iter)|ID),
data=learn_df %>% filter(phase==2),
family="binomial", control = glmerControl(optimizer = "bobyqa"))) ## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(set_size) * scale(stim_iter) + (scale(set_size) *
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 2)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 32506.6 32627.4 -16239.3 32478.6 41236
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.0386 0.1654 0.3001 0.4631 1.3998
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.61266 0.7827
## scale(set_size) 0.11233 0.3352 0.57
## scale(stim_iter) 0.02303 0.1517 0.56 0.36
## scale(set_size):scale(stim_iter) 0.01768 0.1330 0.26 0.24 -0.56
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.02297 0.05113 39.568 <2e-16 ***
## scale(set_size) -0.28678 0.02877 -9.969 <2e-16 ***
## scale(stim_iter) 0.69938 0.02066 33.849 <2e-16 ***
## scale(set_size):scale(stim_iter) -0.25844 0.02082 -12.410 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) 0.298
## scl(stm_tr) 0.410 -0.049
## scl(s_):(_) 0.003 0.436 -0.291# Singular
# summary(pcor_before_after <-
# glmer(correct ~ scale(set_size)*scale(before_after_break) +
# (scale(set_size)*scale(before_after_break)|ID),
# data=end_p1_begin_p2_df,
# family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(pcor_before_after <-
glmer(correct ~ scale(set_size)*scale(before_after_break) +
(scale(set_size) + scale(before_after_break)|ID),
data=end_p1_begin_p2_df,
family="binomial", control = glmerControl(optimizer = "bobyqa"))) ## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ scale(set_size) * scale(before_after_break) + (scale(set_size) +
## scale(before_after_break) | ID)
## Data: end_p1_begin_p2_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 16446.4 16523.6 -8213.2 16426.4 16490
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7814 0.1644 0.3936 0.5720 1.5714
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.421477 0.64921
## scale(set_size) 0.068597 0.26191 0.31
## scale(before_after_break) 0.002658 0.05156 0.34 -0.44
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.41359 0.04558 31.02 <2e-16
## scale(set_size) -0.29171 0.02908 -10.03 <2e-16
## scale(before_after_break) -0.52282 0.02316 -22.58 <2e-16
## scale(set_size):scale(before_after_break) 0.37378 0.02349 15.91 <2e-16
##
## (Intercept) ***
## scale(set_size) ***
## scale(before_after_break) ***
## scale(set_size):scale(before_after_break) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(_) sc(__)
## scal(st_sz) 0.025
## scl(bfr_f_) -0.146 0.229
## scl(_):(__) 0.156 -0.362 -0.277car::vif(pcor_before_after)## scale(set_size)
## 1.174968
## scale(before_after_break)
## 1.105577
## scale(set_size):scale(before_after_break)
## 1.205734car::vif(pcor_phase_2_v2)## scale(set_size) scale(stim_iter)
## 1.245388 1.101492
## scale(set_size):scale(stim_iter)
## 1.357267car::vif(pcor_phase_1_v2)## scale(set_size) scale(stim_iter)
## 1.533658 1.216770
## scale(set_size):scale(stim_iter)
## 1.692973Logistic model predictions are a bit mis-specified but decent
sjPlot::plot_model(pcor_phase_1_v2,
type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") +
ylab("Proportion correct") + ggtitle("Phase 1 regression model predictions")
sjPlot::tab_model(pcor_phase_1_v2)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.71 | 2.55 – 2.88 | <0.001 |
| set size | 0.58 | 0.56 – 0.61 | <0.001 |
| stim iter | 2.92 | 2.78 – 3.06 | <0.001 |
| set size × stim iter | 0.69 | 0.66 – 0.71 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.20 | ||
| τ11 ID.scale(set_size) | 0.05 | ||
| τ11 ID.scale(stim_iter) | 0.10 | ||
| τ11 ID.scale(set_size):scale(stim_iter) | 0.02 | ||
| ρ01 | -0.08 | ||
| 0.90 | |||
| -0.16 | |||
| ICC | 0.10 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.301 / 0.371 | ||
sjPlot::plot_model(pcor_phase_2_v2,
type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") +
ylab("Proportion correct") + ggtitle("Phase 2 regression model predictions")
sjPlot::tab_model(pcor_phase_2_v2)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 7.56 | 6.84 – 8.36 | <0.001 |
| set size | 0.75 | 0.71 – 0.79 | <0.001 |
| stim iter | 2.01 | 1.93 – 2.10 | <0.001 |
| set size × stim iter | 0.77 | 0.74 – 0.80 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.61 | ||
| τ11 ID.scale(set_size) | 0.11 | ||
| τ11 ID.scale(stim_iter) | 0.02 | ||
| τ11 ID.scale(set_size):scale(stim_iter) | 0.02 | ||
| ρ01 | 0.57 | ||
| 0.56 | |||
| 0.26 | |||
| ICC | 0.19 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.136 / 0.299 | ||
pcor_before_after_p <- sjPlot::plot_model(pcor_before_after,
type = "pred", terms = c("before_after_break", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") +
ylab("Proportion correct") + ggtitle("Before to after break regression model predictions")
pcor_before_after_p 
sjPlot::tab_model(pcor_before_after)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 4.11 | 3.76 – 4.49 | <0.001 |
| set size | 0.75 | 0.71 – 0.79 | <0.001 |
| before after break | 0.59 | 0.57 – 0.62 | <0.001 |
|
set size × before after break |
1.45 | 1.39 – 1.52 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.42 | ||
| τ11 ID.scale(set_size) | 0.07 | ||
| τ11 ID.scale(before_after_break) | 0.00 | ||
| ρ01 | 0.31 | ||
| 0.34 | |||
| ICC | 0.13 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.116 / 0.231 | ||
#ggsave("../paper/figs/supp-figs/pcor_before_after_p.png", pcor_before_after_p, height = 5, width=11, dpi=300)Reaction times during the learning phase
Trial 1 so there’s def no S-A-O load on WM
rt_trial_1 <- data.frame(learn_df %>% filter(stim_iter==1) %>%
group_by(set_size, ID) %>% summarize(m=mean(rt)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.rt_trial_1_m <- data.frame(learn_df %>% filter(stim_iter==1) %>%
group_by(set_size) %>% summarize(m=mean(rt)))
rt_tr1_err <- Rmisc::summarySEwithin(rt_trial_1,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_sizeassert(all(round(rt_tr1_err$m, 4)==round(rt_trial_1_m$m, 4)))emp_rt_tr_1 <- ggplot(rt_tr1_err ,
aes(x=set_size, y=m, group=as.factor(set_size), fill=as.factor(set_size))) +
geom_errorbar(aes(x=set_size, ymin=m-se, ymax=m+se), size=1.5, width=.25) +
geom_point(size=6, pch=21, color="black", alpha=.7) +
ga + ap + lp + xlab("Set size") + ylab("") +
ggtitle("At trial 1") + tp + tol#+ ylim(600, 900)rt_ss_si6 <- data.frame(learn_df %>% filter(stim_iter==6) %>%
group_by(set_size, ID) %>% summarize(m=mean(rt)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.rt_ss_si6_m <- data.frame(learn_df %>% filter(stim_iter==6) %>%
group_by(set_size) %>% summarize(m=mean(rt)))
rt_ss_si6_err <- Rmisc::summarySEwithin(rt_ss_si6,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_sizeassert(all(round(rt_ss_si6_err$m, 4)==round(rt_ss_si6_m$m, 4)))emp_rt_si_6 <- ggplot(rt_ss_si6_err ,
aes(x=set_size, y=m, group=as.factor(set_size), fill=as.factor(set_size))) +
geom_errorbar(aes(x=set_size, ymin=m-se, ymax=m+se), size=1.5, width=.25) +
geom_point(size=6, pch=21, color="black", alpha=.7) +
ga + ap + lp + xlab("Set size") + ylab("") +
ggtitle("At stimulus iteration 6") + tp + tolrt_ss_si6_df <- data.frame(learn_df %>% filter(stim_iter==6))rt_ss_si6_df$set_size_indicator <- if_else(rt_ss_si6_df$set_size==1, 1, -1)
# Reduced because model without RE did not converge
summary(rt_ss_si6_df_model <-
lmer(rt ~ set_size_indicator + (1|ID), data=rt_ss_si6_df))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ set_size_indicator + (1 | ID)
## Data: rt_ss_si6_df
##
## REML criterion at convergence: 111694.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5966 -0.6484 -0.0592 0.6007 3.7289
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 7416 86.11
## Residual 41809 204.47
## Number of obs: 8250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 737.512 6.880 593.431 107.20 <2e-16 ***
## set_size_indicator 49.475 4.512 7974.000 10.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## st_sz_ndctr 0.568rt_ss <- data.frame(learn_df %>% group_by(set_size, stim_iter, ID) %>% summarize(m=mean(rt)))## `summarise()` has grouped output by 'set_size', 'stim_iter'. You can override
## using the `.groups` argument.rt_ss_m <- data.frame(learn_df %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(rt))) ## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.rt_ss_err <- Rmisc::summarySEwithin(rt_ss,
measurevar = "m",
withinvars = c("set_size", "stim_iter"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_size, stim_iterassert(all(round(rt_ss_err$m, 4)==round(rt_ss_m$m, 4)))
emp_rt <-
ggplot(rt_ss_err, aes(x=stim_iter, y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_vline(xintercept=c(3, 7, 9), linetype="dotted") +
geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=300, ymax=850, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Reaction time") + tol #+
#ggtitle("") + tp #+ ylim(300, 900)rt_ss_d0_r1 <-
data.frame(learn_df %>% filter(delay==0 & past_reward == 1) %>%
group_by(set_size, stim_iter, ID) %>% summarize(m=mean(rt)))## `summarise()` has grouped output by 'set_size', 'stim_iter'. You can override
## using the `.groups` argument.rt_ss_d0_r1_m <- data.frame(learn_df %>% filter(delay==0 & past_reward == 1) %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(rt)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.rt_ss_err_d0_r1 <- Rmisc::summarySEwithin(rt_ss_d0_r1,
measurevar = "m",
withinvars = c("set_size", "stim_iter"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_size, stim_iterassert(all(round(rt_ss_err_d0_r1$m, 4)==round(rt_ss_err_d0_r1$m, 4)))
emp_rt_d0_r1 <- ggplot(rt_ss_err_d0_r1, aes(x=stim_iter, y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_vline(xintercept=c(2, 5, 7), linetype="dotted") +
geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=5, xmax=8.5, ymin=300, ymax=630, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("") +
ggtitle("When prior reward and no delay") + tp + tol #+ ylim(300, 900) rt_nodrew <- learn_df %>% filter(delay==0 & past_reward == 1)
summary(rt_nodrew_model <-
lmer(rt ~ scale(set_size) + (scale(set_size)|ID), data=rt_nodrew))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) + (scale(set_size) | ID)
## Data: rt_nodrew
##
## REML criterion at convergence: 239851.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1271 -0.6197 -0.1287 0.4544 6.7795
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 4851.9 69.66
## scale(set_size) 480.4 21.92 -0.08
## Residual 20929.7 144.67
## Number of obs: 18678, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 480.237 4.334 274.098 110.81 <2e-16 ***
## scale(set_size) 67.556 1.699 271.451 39.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scal(st_sz) -0.059all_rt <-
emp_rt + emp_rt_d0_r1 + emp_rt_tr_1 + emp_rt_si_6 + plot_annotation(title="Empirical reaction time", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#, size=10)all_rt
#ggsave("../paper/figs/pieces/fig2_emp_rt.png", all_rt, height = 6, width=11, dpi=300)Confirms strong set size effect even at first stim iter — suggesting different proactive strategies
summary(rt_si_1 <-
lmer(rt ~ scale(set_size) + (scale(set_size)|ID), data=learn_df %>% filter(stim_iter==1)))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) + (scale(set_size) | ID)
## Data: learn_df %>% filter(stim_iter == 1)
##
## REML criterion at convergence: 112275.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1250 -0.6697 -0.0766 0.5833 4.1575
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 9725 98.62
## scale(set_size) 1063 32.61 0.31
## Residual 43818 209.33
## Number of obs: 8250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 716.300 6.378 274.026 112.3 <2e-16 ***
## scale(set_size) 37.568 3.030 273.996 12.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scal(st_sz) 0.188sjPlot::tab_model(rt_si_1)| rt | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 716.30 | 703.80 – 728.80 | <0.001 |
| set size | 37.57 | 31.63 – 43.51 | <0.001 |
| Random Effects | |||
| σ2 | 43818.07 | ||
| τ00 ID | 9725.25 | ||
| τ11 ID.scale(set_size) | 1063.36 | ||
| ρ01 ID | 0.31 | ||
| ICC | 0.20 | ||
| N ID | 275 | ||
| Observations | 8250 | ||
| Marginal R2 / Conditional R2 | 0.025 / 0.218 | ||
Regression models in phase 1 and phase 2 show main effects of set size and stim iter and an interaction in both phases
(REs complexity reduce due to convergence failure with full)
summary(rt_p1 <-
lmer(rt ~ scale(set_size)*scale(stim_iter) +
(scale(set_size) + scale(stim_iter)|ID),
data=learn_df %>% filter(phase==1)))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) * scale(stim_iter) + (scale(set_size) +
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 1)
##
## REML criterion at convergence: 558890.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0355 -0.6793 -0.1269 0.5603 4.7845
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 6923.2 83.21
## scale(set_size) 951.4 30.84 0.33
## scale(stim_iter) 428.5 20.70 -0.20 -0.18
## Residual 43257.6 207.98
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 648.359 5.121 274.006 126.61
## scale(set_size) 88.066 2.123 274.000 41.48
## scale(stim_iter) -38.968 1.615 274.014 -24.14
## scale(set_size):scale(stim_iter) 18.135 1.024 40424.000 17.71
## Pr(>|t|)
## (Intercept) <2e-16 ***
## scale(set_size) <2e-16 ***
## scale(stim_iter) <2e-16 ***
## scale(set_size):scale(stim_iter) <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) 0.286
## scl(stm_tr) -0.154 -0.124
## scl(s_):(_) 0.000 0.000 0.000sjPlot::tab_model(rt_p1)| rt | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 648.36 | 638.32 – 658.40 | <0.001 |
| set size | 88.07 | 83.90 – 92.23 | <0.001 |
| stim iter | -38.97 | -42.13 – -35.80 | <0.001 |
| set size × stim iter | 18.13 | 16.13 – 20.14 | <0.001 |
| Random Effects | |||
| σ2 | 43257.63 | ||
| τ00 ID | 6923.16 | ||
| τ11 ID.scale(set_size) | 951.39 | ||
| τ11 ID.scale(stim_iter) | 428.50 | ||
| ρ01 | 0.33 | ||
| -0.20 | |||
| ICC | 0.16 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.157 / 0.293 | ||
sjPlot::plot_model(rt_p1, type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") + ylab("Reaction time") + ggtitle("Phase 1 regression model predictions")
summary(rt_p2 <-lmer(rt ~ scale(set_size)*scale(stim_iter) + (scale(set_size) + scale(stim_iter)|ID), data=learn_df %>% filter(phase==2)))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) * scale(stim_iter) + (scale(set_size) +
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 2)
##
## REML criterion at convergence: 545785.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5886 -0.6438 -0.1083 0.5057 7.5230
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 6462.4 80.39
## scale(set_size) 626.8 25.04 0.28
## scale(stim_iter) 143.2 11.97 -0.12 -0.19
## Residual 31520.5 177.54
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 580.8416 4.9258 273.9913 117.92
## scale(set_size) 57.9437 1.7445 274.0246 33.22
## scale(stim_iter) -40.3142 1.1336 273.9947 -35.56
## scale(set_size):scale(stim_iter) 19.6083 0.8742 40423.9930 22.43
## Pr(>|t|)
## (Intercept) <2e-16 ***
## scale(set_size) <2e-16 ***
## scale(stim_iter) <2e-16 ***
## scale(set_size):scale(stim_iter) <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) 0.234
## scl(stm_tr) -0.075 -0.104
## scl(s_):(_) 0.000 0.000 0.000car::vif(rt_p1)## scale(set_size) scale(stim_iter)
## 1.015515 1.015515
## scale(set_size):scale(stim_iter)
## 1.000000car::vif(rt_p2)## scale(set_size) scale(stim_iter)
## 1.010985 1.010985
## scale(set_size):scale(stim_iter)
## 1.000000sjPlot::tab_model(rt_p2)| rt | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 580.84 | 571.19 – 590.50 | <0.001 |
| set size | 57.94 | 54.52 – 61.36 | <0.001 |
| stim iter | -40.31 | -42.54 – -38.09 | <0.001 |
| set size × stim iter | 19.61 | 17.89 – 21.32 | <0.001 |
| Random Effects | |||
| σ2 | 31520.54 | ||
| τ00 ID | 6462.37 | ||
| τ11 ID.scale(set_size) | 626.77 | ||
| τ11 ID.scale(stim_iter) | 143.23 | ||
| ρ01 | 0.28 | ||
| -0.12 | |||
| ICC | 0.19 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.122 / 0.286 | ||
sjPlot::plot_model(rt_p2, type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") + ylab("Reaction time") + ggtitle("Phase 2 regression model predictions")
Correlation between learning and test phase
learn_ms <- data.frame(learn_df %>%
group_by(ID, set_size) %>% summarize(learn_m=mean(correct)))## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.test_ms <- data.frame(test_df %>%
group_by(ID, set_size) %>% summarize(test_m=mean(correct)))## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.assert(all(test_ms$ID==learn_ms$ID))test_ms$learn_m <- learn_ms$learn_memp_corrs <- ggplot(test_ms,
aes(x=learn_m, y=test_m, color=as.factor(set_size))) +
geom_jitter(alpha=.4, size=2, width=.02) +
ggpubr::stat_cor(method="spearman", size=5, label.y=1.12) + ga + ft + ap +
xlab("Learning") + ylab("Test") + facet_wrap(~ set_size, nrow=1) + tol +
xlim(0, 1) + ylim(0, 1.17) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggtitle("Proportion correct correlation \n Empirical") + tpTake one simulation iter to match number used empirically and not “smooth” by more sims
learn_ms_s <- data.frame(m35_s_learn %>% filter(iter==1) %>% group_by(ID, set_size) %>% summarize(learn_m=mean(corrects)))## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.test_ms_s <- data.frame(m35_s_test %>% filter(iter==1) %>% group_by(ID, set_size) %>% summarize(test_m=mean(corrects)))## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.assert(all(test_ms_s$ID==learn_ms_s$ID))
test_ms_s$learn_m <- learn_ms_s$learn_m
sim_corrs <- ggplot(test_ms_s, aes(x=learn_m, y=test_m, color=as.factor(set_size))) + geom_jitter(alpha=.4, size=2, width=.0) +
ggpubr::stat_cor(method="spearman", size=5, label.y=1.12) + ga + ft + ap +
xlab("Learning") + ylab("Test") + facet_wrap(~ set_size, nrow=1) + tol +
xlim(0, 1) + ylim(0, 1.17) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggtitle("Simulated") + tpA few missing due to jitter and divide by 0
all_corrs <- emp_corrs / sim_corrsall_corrs## Warning: Removed 17 rows containing missing values (`geom_point()`).
#ggsave("../paper/figs/pieces/fig2_perf_corrs.png", all_corrs, height = 6, width=11, dpi=300)Test correct (ie. reward)
Empirical
p1_test_trials <- test_df %>% filter(phase==1)
p2_test_trials <- test_df %>% filter(phase==2)
pcor_p1_test <-
data.frame(p1_test_trials %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_p2_test <-
data.frame(p2_test_trials %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_si6 <- learn_df %>% filter(stim_iter==6) %>%
group_by(set_size, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_p1_test_err <- Rmisc::summarySEwithin(pcor_p1_test,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_sizepcor_si6_err <- Rmisc::summarySEwithin(pcor_si6,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_sizepcor_p2_test_err <- Rmisc::summarySEwithin(pcor_p2_test,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_size# Continuing to sanity check summarSEwithin
pcor_p1_test_m <-
data.frame(p1_test_trials %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p2_test_m <-
data.frame(p2_test_trials %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_si6_m <- learn_df %>% filter(stim_iter==6) %>%
group_by(set_size) %>% summarize(m=mean(correct))
assert(all(round(pcor_p1_test_err$m, 4)==round(pcor_p1_test_m$m, 4)))
assert(all(round(pcor_p2_test_err$m, 4)==round(pcor_p2_test_m$m, 4)))
assert(all(round(pcor_si6_err$m, 4)==round(pcor_si6_m$m, 4)))
emp_p1_test <- ggplot(pcor_p1_test_err, aes(x=set_size, y=m, fill=set_size)) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=m-se, ymax=m+se), width=.2) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(0, 1) +
tp + ggtitle("Test phase 1")
emp_si6 <- ggplot(pcor_si6_err, aes(x=set_size, y=m, fill=set_size)) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=m-se, ymax=m+se), width=.2) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) +
tp + ggtitle("Stimulus iteration 6")
emp_p2_test <- ggplot(pcor_p2_test_err, aes(x=set_size, y=m, fill=set_size)) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=m-se, ymax=m+se), width=.2) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) + tp +
ggtitle("Test phase 2")U-shaped pattern evident at test phase 1, which replicates at the beginning of the next learning phase (stimulus iteration 6) before there has been a chance for further learning (thus making it effectively another test trial)
The effect then appears to reduce substantially by test phase 2 — with stim iter 1 still decreased, but set size 2 if anything higher than the others
emps_u_plot <-
emp_p1_test + emp_si6 + emp_p2_test + plot_annotation(title="Empirical", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#,emps_u_plot
Simulation plots
pcor_p1_test_sim_m35 <-
data.frame(m35_s_test %>% filter(phase==1) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p1_test_sim_m35_iters <-
data.frame(m35_s_test %>% filter(phase==1) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_si6_sim_m35 <-
data.frame(m35_s_learn %>% filter(stim_iter==6) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_si6_sim_m35_iters <-
data.frame(m35_s_learn %>% filter(stim_iter==6) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_p2_test_sim_m35 <-
data.frame(m35_s_test %>% filter(phase==2) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p2_test_sim_m35_iters <-
data.frame(m35_s_test %>% filter(phase==2) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.sim_p1_test_m35 <-
ggplot(pcor_p1_test_sim_m35, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p1_test_sim_m35_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(0, 1) +
tp + ggtitle("Test phase 1")
sim_si6_m35 <- ggplot(pcor_si6_sim_m35, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_si6_sim_m35_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) +
tp + ggtitle("Stimulus iteration 6")
sim_p2_test_m35 <- ggplot(pcor_p2_test_sim_m35, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p2_test_sim_m35_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) + tp +
tp + ggtitle("Test phase 2")sims_u_plot <-
sim_p1_test_m35 + sim_si6_m35 + sim_p2_test_m35 + plot_annotation(title="Simulated", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#,sims_u_plot
# ggsave("../paper/figs/pieces/fig3_empU.png", emps_u_plot, height = 4, width=11, dpi=300)
#ggsave("../paper/figs/pieces/fig3_simU.png", sims_u_plot, height = 4, width=11, dpi=300)summary(pcor_p1_test_quad <-
glmer(correct ~ poly(set_size, 2) +
(poly(set_size, 2)|ID),
data=test_df %>% filter(phase==1),
family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ poly(set_size, 2) + (poly(set_size, 2) | ID)
## Data: test_df %>% filter(phase == 1)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 18853.1 18922.5 -9417.5 18835.1 16491
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9057 -0.8769 0.4293 0.6433 1.9108
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.6646 0.8152
## poly(set_size, 2)1 1999.3502 44.7141 0.15
## poly(set_size, 2)2 2394.8985 48.9377 -0.13 -0.12
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95321 0.05283 18.043 < 2e-16 ***
## poly(set_size, 2)1 13.42438 3.24378 4.138 3.50e-05 ***
## poly(set_size, 2)2 -26.23676 3.43764 -7.632 2.31e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2)1
## ply(st_,2)1 0.096
## ply(st_,2)2 -0.099 -0.065Predictions look reasonable and do fit a clear inverted U
sjPlot::tab_model(pcor_p1_test_quad)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.59 | 2.34 – 2.88 | <0.001 |
| set size [1st degree] | 676291.41 | 1172.25 – 390163629.22 | <0.001 |
| set size [2nd degree] | 0.00 | 0.00 – 0.00 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.66 | ||
| τ11 ID.poly(set_size, 2)1 | 1999.35 | ||
| τ11 ID.poly(set_size, 2)2 | 2394.90 | ||
| ρ01 | 0.15 | ||
| -0.13 | |||
| ICC | 0.22 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.012 / 0.230 | ||
pcor_p1_test_quad_p <- sjPlot::plot_model(pcor_p1_test_quad, type = "pred",
terms = c("set_size [all]")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("Proportion correct") + ggtitle("Quadratic model regression predictions \n Phase 1 test")
pcor_p1_test_quad_p
Stimulus iteration 6 replicates effect
summary(pcor_si6_quad <-
glmer(correct ~ poly(set_size, 2) +
(poly(set_size, 2)|ID),
data=learn_df %>% filter(stim_iter==6),
family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ poly(set_size, 2) + (poly(set_size, 2) | ID)
## Data: learn_df %>% filter(stim_iter == 6)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 9777.1 9840.3 -4879.6 9759.1 8241
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9925 -1.0132 0.4758 0.6605 1.4247
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4794 0.6924
## poly(set_size, 2)1 414.1219 20.3500 0.28
## poly(set_size, 2)2 366.8434 19.1532 -0.48 0.62
## Number of obs: 8250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.91428 0.04948 18.478 < 2e-16 ***
## poly(set_size, 2)1 6.44043 2.43970 2.640 0.00829 **
## poly(set_size, 2)2 -22.80614 2.48653 -9.172 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2)1
## ply(st_,2)1 0.122
## ply(st_,2)2 -0.237 0.204pcor_si6_quad_p <- sjPlot::plot_model(pcor_si6_quad,
type = "pred", terms = c("set_size")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("Proportion correct") +
ggtitle("Quadratic model regression predictions \n Stim iter 6")## Model contains splines or polynomial terms. Consider using
## `terms="set_size [all]"` to get smooth plots. See also package-vignette
## 'Marginal Effects at Specific Values'.pcor_si6_quad_p 
sjPlot::tab_model(pcor_si6_quad)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.49 | 2.26 – 2.75 | <0.001 |
| set size [1st degree] | 626.67 | 5.25 – 74768.23 | 0.008 |
| set size [2nd degree] | 0.00 | 0.00 – 0.00 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.48 | ||
| τ11 ID.poly(set_size, 2)1 | 414.12 | ||
| τ11 ID.poly(set_size, 2)2 | 366.84 | ||
| ρ01 | 0.28 | ||
| -0.48 | |||
| ICC | 0.15 | ||
| N ID | 275 | ||
| Observations | 8250 | ||
| Marginal R2 / Conditional R2 | 0.017 / 0.163 | ||
Neither linear nor quadratic signif effects
summary(pcor_p2_test_quad <-
glmer(correct ~ poly(set_size, 2) +
(poly(set_size, 2)|ID),
data=test_df %>% filter(phase==2),
family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ poly(set_size, 2) + (poly(set_size, 2) | ID)
## Data: test_df %>% filter(phase == 2)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 12074.6 12144.0 -6028.3 12056.6 16491
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4912 0.1678 0.2476 0.3975 1.4058
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.513 1.23
## poly(set_size, 2)1 3100.386 55.68 0.02
## poly(set_size, 2)2 2531.491 50.31 0.06 -0.15
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3021 0.0814 28.280 <2e-16 ***
## poly(set_size, 2)1 -3.1372 3.4343 -0.913 0.361
## poly(set_size, 2)2 -4.5952 3.7001 -1.242 0.214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2)1
## ply(st_,2)1 0.010
## ply(st_,2)2 0.025 0.069sjPlot::tab_model(pcor_p2_test_quad)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 10.00 | 8.52 – 11.72 | <0.001 |
| set size [1st degree] | 0.04 | 0.00 – 36.38 | 0.361 |
| set size [2nd degree] | 0.01 | 0.00 – 14.25 | 0.214 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.51 | ||
| τ11 ID.poly(set_size, 2)1 | 3100.39 | ||
| τ11 ID.poly(set_size, 2)2 | 2531.49 | ||
| ρ01 | 0.02 | ||
| 0.06 | |||
| ICC | 0.36 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.361 | ||
pcor_p2_quad_p <- sjPlot::plot_model(pcor_p2_test_quad,
type = "pred", terms = c("set_size [all]")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("Proportion correct") +
ggtitle("Quadratic model regression \npredictions \n Phase 2 test")
pcor_p2_quad_p
pcor_p1_test_quad_p + pcor_si6_quad_p +pcor_p2_quad_p 
pcor_p2_quad_p
Using the RL off parameter to capture individual differences in blunting underlying the inverted U
rl_off_p <-
ggplot(m35, aes(x=rl_off)) +
geom_vline(xintercept=median(m35$rl_off), size=4) +
geom_histogram(fill="white", color="black") + ga + ap + ylab("") +
ggtitle(TeX('$\\RL^{off}')) + tp + xlab("")rl_off_p## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#ggsave("../paper/figs/pieces/rl_off.png", rl_off_p, height = 4, width=11, dpi=300)Dividing up empirically into high vs. low blunters
little_coop <- data.frame(m35 %>% filter(rl_off < median(m35$rl_off)))$ID
high_coop <- data.frame(m35 %>% filter(rl_off > median(m35$rl_off)))$IDp1_test_emp_trials_m35_lc <- test_df %>% filter(phase==1 & ID %in% little_coop)
pcor_p1_test_emp_m35_lc <-
data.frame(p1_test_emp_trials_m35_lc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p1_test_emp_m35_lc_ID <-
data.frame(p1_test_emp_trials_m35_lc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.p2_test_emp_trials_m35_lc <- test_df %>% filter(phase==2 & ID %in% little_coop)
pcor_p2_test_emp_m35_lc <-
data.frame(p2_test_emp_trials_m35_lc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p2_test_emp_m35_lc_ID <-
data.frame(p2_test_emp_trials_m35_lc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.si6_emp_trials_m35_lc <- learn_df %>% filter(stim_iter==6 & ID %in% little_coop)
pcor_si6_m35_lc <-
data.frame(si6_emp_trials_m35_lc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_si6_m35_lc_ID <-
data.frame(si6_emp_trials_m35_lc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.emp_p1_test_m35_lc <- ggplot(pcor_p1_test_emp_m35_lc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p1_test_emp_m35_lc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) +
tp + ggtitle("Test phase 1")
emp_p2_test_m35_lc <- ggplot(pcor_p2_test_emp_m35_lc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p2_test_emp_m35_lc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) +
tp + ggtitle("Test phase 2")
emp_si6_lc <-
ggplot(pcor_si6_m35_lc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_si6_m35_lc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) +
tp + ggtitle("Stimulus iteration 6")low_blunt <- emp_p1_test_m35_lc + emp_si6_lc + emp_p2_test_m35_lcHigh blunters
p1_test_emp_trials_m35_hc <- test_df %>% filter(phase==1 & ID %in% high_coop)
pcor_p1_test_emp_m35_hc <-
data.frame(p1_test_emp_trials_m35_hc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p1_test_emp_m35_hc_ID <-
data.frame(p1_test_emp_trials_m35_hc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.p2_test_emp_trials_m35_hc <- test_df %>% filter(phase==2 & ID %in% high_coop)
pcor_p2_test_emp_m35_hc <-
data.frame(p2_test_emp_trials_m35_hc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p2_test_emp_m35_hc_ID <-
data.frame(p2_test_emp_trials_m35_hc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.si6_emp_trials_m35_hc <- learn_df %>% filter(stim_iter==6 & ID %in% high_coop)
pcor_si6_m35_hc <-
data.frame(si6_emp_trials_m35_hc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_si6_m35_hc_ID <-
data.frame(si6_emp_trials_m35_hc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.emp_p1_test_m35_hc <- ggplot(pcor_p1_test_emp_m35_hc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p1_test_emp_m35_hc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) #+
#tp + ggtitle("Phase 1")
emp_p2_test_m35_hc <- ggplot(pcor_p2_test_emp_m35_hc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p2_test_emp_m35_hc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) #+
# tp + ggtitle("Phase 2")
emp_si6_hc <-
ggplot(pcor_si6_m35_hc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_si6_m35_hc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) #+
#tp + ggtitle("Stimulus iteration 6")high_blunt <- emp_p1_test_m35_hc + emp_si6_hc + emp_p2_test_m35_hclow_vs_high <- low_blunt/high_bluntlow_vs_high
#ggsave("../paper/figs/pieces/fig4_lohiblunt.png", low_vs_high, height = 9, width=11, dpi=300)Put RL off into the test df
unique_ids <- unique(test_df$ID)
for (i in 1:length(unique_ids)) {
test_df[test_df$ID==unique_ids[i], "rl_off"] <- m35[m35$ID==unique_ids[i], "rl_off"]
}
# Spot check
# test_df %>% filter(ID == 24) %>% select(rl_off)
# m35 %>% filter(ID == 24) %>% select(rl_off)
# test_df %>% filter(ID == 85) %>% select(rl_off)
# m35 %>% filter(ID == 85) %>% select(rl_off)
# test_df %>% filter(ID == 2) %>% select(rl_off)
# m35 %>% filter(ID == 2) %>% select(rl_off)
# test_df %>% filter(ID == 219) %>% select(rl_off)
# m35 %>% filter(ID == 219) %>% select(rl_off)m_rlo <- mean(m35$rl_off)
sd_up <- m_rlo +sd(m35$rl_off) # (.5*sd(m35$rl_off))
sd_down <- m_rlo - sd(m35$rl_off)# (.5*sd(m35$rl_off))
cat("\nMean", m_rlo, "\nOne SD up", sd_up, "\nOne SD down", sd_down)##
## Mean 0.544767
## One SD up 0.8678609
## One SD down 0.221673#quantile(m35$rl_off, seq(1/3, 1, 1/3))#ggsave("../paper/figs/pieces/fig4_mod_rl_off_p.png", mod_rl_off_p, height = 4, width=8, dpi=300)summary(pcor_p2_moderation <-
glmer(correct ~ scale(set_size)*scale(rl_off) +
(scale(set_size)*scale(rl_off)|ID),
data=test_df %>% filter(phase==2),
family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(set_size) * scale(rl_off) + (scale(set_size) *
## scale(rl_off) | ID)
## Data: test_df %>% filter(phase == 2)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 12009.4 12117.4 -5990.7 11981.4 16486
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.6425 0.1503 0.2565 0.4166 1.5443
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.20709 1.0987
## scale(set_size) 0.04156 0.2039 0.19
## scale(rl_off) 0.25019 0.5002 -0.34 0.41
## scale(set_size):scale(rl_off) 0.06133 0.2476 0.68 0.41 0.41
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.32737 0.08506 27.362 < 2e-16 ***
## scale(set_size) -0.03490 0.04117 -0.848 0.397
## scale(rl_off) -0.54108 0.09283 -5.829 5.59e-09 ***
## scale(set_size):scale(rl_off) 0.45056 0.04696 9.594 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(s_) scl(r_)
## scal(st_sz) -0.008
## scal(rl_ff) -0.306 0.354
## scl(s_):(_) 0.352 -0.404 -0.038car::vif(pcor_p2_moderation)## scale(set_size) scale(rl_off)
## 1.384028 1.160176
## scale(set_size):scale(rl_off)
## 1.212803mod_rl_off_p <- sjPlot::plot_model(pcor_p2_moderation, type = "pred",
terms = c("set_size [all]", "rl_off [.2217, .5447, .8678]")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("") +
ggtitle("Moderation of set size effect in final test phase \n by RL blunting") +
scale_color_manual(values=c("blue", "gray40", "red"), labels=c("Low", "Medium", "High")) +
scale_fill_manual(values=c("blue", "gray40", "red"), labels=c("Low", "Medium", "High"))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.mod_rl_off_p
sjPlot::tab_model(pcor_p2_moderation)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 10.25 | 8.68 – 12.11 | <0.001 |
| set size | 0.97 | 0.89 – 1.05 | 0.397 |
| rl off | 0.58 | 0.49 – 0.70 | <0.001 |
| set size × rl off | 1.57 | 1.43 – 1.72 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.21 | ||
| τ11 ID.scale(set_size) | 0.04 | ||
| τ11 ID.scale(rl_off) | 0.25 | ||
| τ11 ID.scale(set_size):scale(rl_off) | 0.06 | ||
| ρ01 | 0.19 | ||
| -0.34 | |||
| 0.68 | |||
| ICC | 0.32 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.093 / 0.385 | ||
Low vs. high low set sizes percents reported in paper
mean(pcor_p1_test_emp_m35_hc[1:3, "m"])## [1] 0.5612328mean(pcor_p1_test_emp_m35_lc[1:3, "m"])## [1] 0.729927mean(pcor_si6_m35_hc[1:3, "m"])## [1] 0.5788727mean(pcor_si6_m35_lc[1:3, "m"])## [1] 0.7226277mean(pcor_p2_test_emp_m35_hc[1:3, "m"])## [1] 0.7650041mean(pcor_p2_test_emp_m35_lc[1:3, "m"])## [1] 0.915146Low vs. high high set sizes
mean(pcor_p1_test_emp_m35_hc[4:5, "m"])## [1] 0.6972628mean(pcor_p1_test_emp_m35_lc[4:5, "m"])## [1] 0.7048814mean(pcor_si6_m35_hc[4:5, "m"])## [1] 0.6803832mean(pcor_si6_m35_lc[4:5, "m"])## [1] 0.7076642mean(pcor_p2_test_emp_m35_hc[4:5, "m"])## [1] 0.8436131mean(pcor_p2_test_emp_m35_lc[4:5, "m"])## [1] 0.8611314Checking low vs. high asymptotic performance — set sizes 1 and 2 are very comparable
lc_asym <- learn_df %>% filter(stim_iter==10 & ID %in% little_coop) %>% group_by(set_size) %>% summarize(m=mean(correct))hc_asym <- learn_df %>% filter(stim_iter==10 & ID %in% high_coop) %>% group_by(set_size) %>% summarize(m=mean(correct))mean(unlist(lc_asym[1:3, "m"]))## [1] 0.9578264mean(unlist(hc_asym[1:3, "m"]))## [1] 0.9551906mean(lc_asym[1:3, "m"])## Warning in mean.default(lc_asym[1:3, "m"]): argument is not numeric or logical:
## returning NA## [1] NARL off is not much correlated with WM pars
ComparePars(m35$rl_off, m35$phi, use_identity_line = 0)
ComparePars(m35$rl_off, m35$kappa, use_identity_line = 0)
ComparePars(m35$rl_off, m35$rho, use_identity_line = 0)
Incorrect trials
Neutral (over worst) preference at learning and test
A key change in this task from the classic RLWM was having punishment | neutral | reward bandit arms instead of just reward | neutral | neutral. But it’s possible that pts would focus on maximizing reward and not make much distinguish punishment vs. neutral. The below examine neutral > punishment preference during learning and test (retention).
learn_error_df <- learn_df %>% filter(correct==0)
np_summs_m <- learn_error_df %>% group_by(stim_iter) %>%
summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))
# Within-subject adjusted mean and raw don't agree presumably because of amount of
# variation given sparsity of errors/subject so will just use the raw means for plotting
# so can do accurate compare to model preds
# np_si <- Rmisc::summarySEwithin(np_summs,
# measurevar = "neutral_pref",
# withinvars = c("stim_iter"),
# idvar = "ID")
# These do disagree before this decimal point
# round(np_summs_m$neutral_pref, 4)
# round(np_si$neutral_pref, 4)
emp_np_si_plot <-
ggplot(np_summs_m, aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
#geom_errorbar(aes(ymin=neutral_pref-se, ymax=neutral_pref+se), width=.2) +
geom_bar(data=np_summs_m,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter)),
stat="identity", color="black", alpha=.85) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, fill="gray57") +
annotate("text", x=3, y=.20, label="Phase 1", size=8) +
annotate("text", x=7.5, y=.20, label="Phase 2", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp +
ggtitle("Empirical") + scale_color_gradient2() + ylim(-.1, .23)
emp_np_si_plot
#move_layers(emp_np_si_plot, "GeomPoint", position = "top")Capture by sims
np_si_sims_summs <- m35_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_si_sims_summs_var <- m35_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter, iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))## `summarise()` has grouped output by 'stim_iter'. You can override using the
## `.groups` argument.m35_np_si_plot <- ggplot(np_si_sims_summs,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_jitter(data=np_si_sims_summs_var, aes(x=stim_iter, y=neutral_pref), height=0,
width=.2, pch=21) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, pch=21) +
annotate("text", x=3, y=.20, label="Phase 1", size=8) +
annotate("text", x=7.5, y=.20, label="Phase 2", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp +
ylim(-.05, .23) +
ggtitle("Simulated") + scale_color_gradient2() +
scale_x_continuous(breaks=seq(1, 10, 1)) + ylim(-.1, .23)## Warning in annotate("rect", xmin = 5.5, xmax = 10.5, ymin = 0, ymax = 0.22, :
## Ignoring unknown parameters: `shape`## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.m35_np_si_plot
Simulation w no punish bonus off
m35_luu <- read.csv("../model_res/sims/SIM_RunRLWMPRewLesionNotPun67064.csv") # Updated non-bug 8-1-24 model np_si_sims_summs_luu <- m35_luu %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_si_sims_summs_var_luu <- m35_luu %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))## `summarise()` has grouped output by 'stim_iter'. You can override using the
## `.groups` argument.m35_np_si_luu_plot <-
ggplot(np_si_sims_summs_luu,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_jitter(data=np_si_sims_summs_var_luu, aes(x=stim_iter, y=neutral_pref), pch=21, width=.2) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, fill="gray57") +
annotate("text", x=3, y=.20, label="Phase 1", size=8) +
annotate("text", x=7.5, y=.20, label="Phase 2", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp + ylim(-.1, .23) +
ggtitle(TeX("Sim.: Non-pun. bonus = 0")) + scale_color_gradient2() +
scale_x_continuous(breaks=seq(1, 10, 1)) + theme(plot.title = element_text(size = 18))
#theme(plot.title = element_text(vjust = -2.5))
m35_np_si_luu_plot
np_si <- emp_np_si_plot + m35_np_si_plot + m35_np_si_luu_plotnp_si
#ggsave("../paper/figs/pieces/fig5_neut-pref_emp-and-sims.png", np_si, height = 6, width=12, dpi=300)# Didn't converge, so do Bayesian version
# summary(pneut_phase1_le <- glmer(neutral ~ scale(stim_iter) + (0 + scale(stim_iter) |ID),
# data=learn_error_df %>% filter(phase==1), family="binomial", control = glmerControl(optimizer = "bobyqa")))
# Now runnning in brms_s script
# pneut_time <- brm(
# neutral ~ scale(stim_iter) + (1 |ID),
# data = learn_error_df %>% filter(phase==1),
# family = bernoulli(link = "logit"),
# warmup=2e3,
# iter=4e3,
# chains=5,
# cores=5,
# control= list(adapt_delta = 0.9))Note: brms files not shared in public-data due to file
size but can be reproduced using brms-s.R
np_p1_brms <-
read.csv("../model_res/brms_res/neutral_pref_time_effect_phase1__66703__.csv")
np_p2_brms <-
read.csv("../model_res/brms_res/neutral_pref_time_effect_phase2__28768__.csv")hist(np_p1_brms$b_scalestim_iter, breaks=100)
hist(np_p2_brms$b_scalestim_iter, breaks=100)
Testing effects
test_error_df <- test_df %>% filter(correct==0)
# np_test_summs <- test_error_df %>% group_by(ID, phase) %>%
# summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_test_summs_m <- test_error_df %>% group_by(phase) %>%
summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))
# np_test_si <- Rmisc::summarySEwithin(np_test_summs,
# measurevar = "neutral_pref",
# withinvars = c("phase"),
# idvar = "ID")
#
# # Within-subject adjusted mean and raw don't agree presumably because of amount of
# # variation given sparsity of errors/subject so will just use the raw means for plotting
# # so can do accurate compare to model preds
# np_test_si_sew <- Rmisc::summarySEwithin(np_test_summs,
# measurevar = "neutral_pref",
# withinvars = c("phase"),
# idvar = "ID")
#
# round(np_test_summs_m$neutral_pref, 4)
# round(np_test_si_sew$neutral_pref, 4)In contrast to during learning, participants do not appear to retain much ability to avoid the wost stimuli during test
np_test_summs_m## # A tibble: 2 × 5
## phase mw mn n neutral_pref
## <int> <dbl> <dbl> <int> <dbl>
## 1 1 0.494 0.506 5129 0.0119
## 2 2 0.488 0.512 2458 0.0236emp_np_si_test_plot <-
ggplot(np_test_summs_m, aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
#geom_errorbar(aes(ymin=neutral_pref-se, ymax=neutral_pref+se), width=.2) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(#"Neutral preference during test
"Empirical") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
emp_np_si_test_plot
# ggsave("../paper/figs/pieces/fig6_emp_np_si.png", emp_np_si_test_plot, height = 4, width=5, dpi=300)And statistically there is no significant ability to pick neutral over worst, or better ability to do so by the second versus first phase.
(These are not singular issue so didn’t need to fit Bayesian model)
summary(pneut_test <- glmer(neutral ~ 1 + (1|ID),
data=test_error_df, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: neutral ~ 1 + (1 | ID)
## Data: test_error_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 10474.7 10488.6 -5235.4 10470.7 7585
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4510 -0.9779 0.7725 0.9626 1.3129
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.09459 0.3076
## Number of obs: 7587, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.03924 0.03108 1.263 0.207summary(pneut_test_phase <- glmer(neutral ~ phase + (phase|ID),
data=test_error_df, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: neutral ~ phase + (phase | ID)
## Data: test_error_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 10477.0 10511.6 -5233.5 10467.0 7582
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6292 -0.9788 0.7230 0.9641 1.3537
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.12691 0.3562
## phase 0.06782 0.2604 -0.61
## Number of obs: 7587, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.02154 0.07496 -0.287 0.774
## phase 0.04757 0.05488 0.867 0.386
##
## Correlation of Fixed Effects:
## (Intr)
## phase -0.910Bayesian version for more direct comparison to learning where only Bayesian model converged
np_test_brms <- read.csv("../model_res/brms_res/neutral_pref_test_effect__33330__.csv")hist(np_test_brms$b_Intercept, breaks=100)
However store the frequentist REs for correlating with test diffs
test_res <- data.frame(ranef(pneut_test)$ID)The lack of retained ability to avoid punishment over neutral is unsurprising because the RL learning rate from negative PEs fit much lower than that from positive PEs
cat("\nAlpha for pos. PEs quantiles: \n")##
## Alpha for pos. PEs quantiles:round(quantile(m35$alpha_pos), 6)## 0% 25% 50% 75% 100%
## 0.000767 0.005559 0.007560 0.010339 0.048076cat("\nAlpha for neg. PEs quantiles: \n")##
## Alpha for neg. PEs quantiles:round(quantile(m35$alpha_neg), 6)## 0% 25% 50% 75% 100%
## 0.000000 0.000000 0.000122 0.002330 0.023870Take log given strong negative skew
hist(m35$alpha_neg, breaks=100)
hist(log(m35$alpha_neg), breaks=100)
hist(m35$alpha_pos, breaks=100)
hist(log(m35$alpha_pos), breaks=100)
t.test(log(m35$alpha_neg), log(m35$alpha_pos), paired = TRUE)##
## Paired t-test
##
## data: log(m35$alpha_neg) and log(m35$alpha_pos)
## t = -16.124, df = 274, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -7.853473 -6.144423
## sample estimates:
## mean difference
## -6.998948m35_test_sim_error_df <- m35_s_test %>% filter(corrects==0)
np_test_summs_m35_sim <- m35_test_sim_error_df %>% group_by(phase) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_test_summs_m35_sim_var <- m35_test_sim_error_df %>% group_by(phase, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))## `summarise()` has grouped output by 'phase'. You can override using the
## `.groups` argument.And the model correctly predicts poor retained ability to avoid punishment over neutral overall
sim_m35_np_si_test_plot <-
ggplot(np_test_summs_m35_sim , aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_jitter(data=np_test_summs_m35_sim_var, height=0, width=.2, size=3, alpha=.8) +
geom_bar(stat="identity", color="black", alpha=.8) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(#"Neutral preference during Test
"Simulated") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
sim_m35_np_si_test_plot
Whereas a model that sets \(\alpha_{neg} = \alpha_{pos}\) incorrectly predicts a retained neutral preference at a similar level to what is observed empirically at the highest points during learning
m35_s_lesionalphaneg <-
read.csv("../model_res/sims/SIM_RunRLWMPRewLesionAlphaNeg27102.csv")m35_test_sim_error_df_lan <- m35_s_lesionalphaneg %>% filter(corrects==0)
np_test_summs_m35_sim_lan <- m35_test_sim_error_df_lan %>% group_by(phase) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_test_summs_m35_sim_var_lan <- m35_test_sim_error_df_lan %>% group_by(phase, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))## `summarise()` has grouped output by 'phase'. You can override using the
## `.groups` argument.sim_m35_np_si_test_plot_lan <-
ggplot(np_test_summs_m35_sim_lan , aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_jitter(data=np_test_summs_m35_sim_var_lan, height=0, width=.2, size=3, alpha=.8) +
geom_bar(stat="identity", color="black", alpha=.8) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(TeX("Simulated: $\\alpha^{-} = \\alpha^{+}$")) +
scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
#ggtitle("Neutral preference during Test \nSimulated equating TeX('$\\alpha^{-}')") +
sim_m35_np_si_test_plot_lan
test_np_plots <-
emp_np_si_test_plot + sim_m35_np_si_test_plot #+ sim_m35_np_si_test_plot_lantest_np_plots
ggsave("../paper/figs/pieces/fig6_TEST_neut-pref_emp-and-sims.png", test_np_plots, height = 6, width=12, dpi=300)Code if alpha pos is greater
m35$alpha_pos_greater <- if_else(m35$alpha_pos > m35$alpha_neg, 1, 0)table(m35$alpha_pos_greater)##
## 0 1
## 17 258table(m35$alpha_pos_greater)[2]/sum(table(m35$alpha_pos_greater))## 1
## 0.9381818alpha_comp <- ggplot(m35,
aes(x=alpha_neg, alpha_pos, fill=as.factor(alpha_pos_greater))) +
geom_point(size=4, pch=21) +
ga + ap + tp +
stp + tol + ylim(0, .05) + xlim(0, .05) +
geom_line(aes(y=alpha_neg)) +
xlab(TeX("$\\alpha^{-} ")) +
ylab(TeX("$\\alpha^{+} ")) + scale_fill_manual(values=c("purple", "orange")) +
theme(axis.title = element_text(size=50))
# ggtitle(model_char, subtitle=str) +
# ylab(ychar) + xlab(xchar)
alpha_comp
alpha_and_sim <- alpha_comp + sim_m35_np_si_test_plot_lan# ggsave("../paper/figs/pieces/fig5_alpha_and_sim.png", alpha_and_sim, height = 6, width=12, dpi=300)However, there are some individual differences — with those fitting a higher learning from negative PEs showing better ability to retain a neutral preference (and note this correlation is expected to be attenuated due to measurement error due to the relatively poor parameter recovery for the learning rate for negative PEs shown earlier)
hist(m35$alpha_neg)
ComparePars(m35$alpha_neg, test_res$X.Intercept.,
model_char="", xchar = "Learning rate from negative RPEs",
ychar = "Avoidance of punishment \n over neutral at test",
use_identity_line = 0)
ComparePars(log(m35$alpha_neg), test_res$X.Intercept.,
model_char="", xchar = "Log (learning rate from negative RPEs)",
ychar = "Avoidance of punishment \n over neutral at test",
use_identity_line = 0)
Individual differences - depression/anxiety and rumination
# qdf %>% filter(catch_q_1 != 1) %>% select(ID)
# qdf %>% filter(catch_q_2 != 1) %>% select(ID)
qdf <- read.csv("../data/questionnaire_df.csv")
demogs <- read.csv("../data/demogs_deident.csv")
if (all(demogs$ID==qdf$ID)) qdf$group <- demogs$group
qdf <- qdf %>% filter(!ID==255)
demogs <- demogs %>% filter(!ID==255)
learn_df <- learn_df %>% filter(!ID==255)
test_df <- test_df %>% filter(!ID==255)Tiny amount of missing cases so fine to mean impute
length(which(is.na(qdf %>% select(contains("BDI"))))) # 3 missing BDI ## [1] 3length(which(is.na(qdf %>% select(contains("GAD"))))) # 4 GAD ## [1] 4length(which(is.na(qdf %>% select(contains("RRS"))))) # 3 RRS## [1] 3qdf_before_imp <- qdffor (i in 1:nrow(qdf)) {
if (any(is.na(qdf[i, grep("BDI", names(qdf))]))) {
this_row <- qdf[i, grep("BDI", names(qdf))]
na_idx <- which(is.na(qdf[i, grep("BDI", names(qdf))]))
this_row[na_idx] <- mean(as.numeric(this_row), na.rm=TRUE)
qdf[i, grep("BDI", names(qdf))] <- this_row
}
}
for (i in 1:nrow(qdf)) {
if (any(is.na(qdf[i, grep("GAD", names(qdf))]))) {
this_row <- qdf[i, grep("GAD", names(qdf))]
na_idx <- which(is.na(qdf[i, grep("GAD", names(qdf))]))
this_row[na_idx] <- mean(as.numeric(this_row), na.rm=TRUE)
qdf[i, grep("GAD", names(qdf))] <- this_row
}
}
for (i in 1:nrow(qdf)) {
if (any(is.na(qdf[i, grep("RRS", names(qdf))]))) {
this_row <- qdf[i, grep("RRS", names(qdf))]
na_idx <- which(is.na(qdf[i, grep("RRS", names(qdf))]))
this_row[na_idx] <- mean(as.numeric(this_row), na.rm=TRUE)
qdf[i, grep("RRS", names(qdf))] <- this_row
}
}Given tiny amount of missing cases so shouldn’t skew values much, doing NA rm for pres
length(which(is.na(qdf %>% select(contains("BDI"))))) # 3 missing BDI ## [1] 0length(which(is.na(qdf %>% select(contains("GAD"))))) # 4 GAD ## [1] 0length(which(is.na(qdf %>% select(contains("RRS"))))) # 3 RRS## [1] 0Sanity check unchanged
# ComparePars(rowSums(qdf[, grep("RRS", names(qdf))]),
# rowSums(qdf_before_imp[, grep("RRS", names(qdf_before_imp))]))
# ComparePars(rowSums(qdf[, grep("BDI", names(qdf))]),
# rowSums(qdf_before_imp[, grep("BDI", names(qdf_before_imp))]))
# ComparePars(rowSums(qdf[, grep("GAD", names(qdf))]),
# rowSums(qdf_before_imp[, grep("GAD", names(qdf_before_imp))]))Reliability stats
psych::alpha(qdf[, grep("BDI", names(qdf))])##
## Reliability analysis
## Call: psych::alpha(x = qdf[, grep("BDI", names(qdf))])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.96 0.96 0.97 0.54 25 0.0033 0.68 0.66 0.54
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.95 0.96 0.97
## Duhachek 0.95 0.96 0.97
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## BDI.1 0.96 0.96 0.97 0.54 24 0.0034 0.0107 0.54
## BDI.2 0.96 0.96 0.97 0.54 23 0.0035 0.0102 0.54
## BDI.3 0.96 0.96 0.97 0.54 23 0.0036 0.0099 0.54
## BDI.4 0.96 0.96 0.97 0.54 23 0.0035 0.0102 0.54
## BDI.5 0.96 0.96 0.97 0.55 24 0.0034 0.0111 0.55
## BDI.6 0.96 0.96 0.97 0.55 24 0.0034 0.0107 0.55
## BDI.7 0.96 0.96 0.97 0.54 23 0.0036 0.0099 0.54
## BDI.8 0.96 0.96 0.97 0.54 23 0.0036 0.0101 0.54
## BDI.9 0.96 0.96 0.97 0.55 24 0.0034 0.0103 0.54
## BDI.10 0.96 0.96 0.97 0.55 24 0.0034 0.0107 0.56
## BDI.11 0.96 0.96 0.97 0.54 24 0.0034 0.0111 0.55
## BDI.12 0.96 0.96 0.97 0.54 23 0.0036 0.0100 0.54
## BDI.13 0.96 0.96 0.97 0.54 24 0.0035 0.0109 0.54
## BDI.14 0.96 0.96 0.97 0.54 23 0.0036 0.0098 0.54
## BDI.15 0.96 0.96 0.97 0.54 23 0.0035 0.0104 0.54
## BDI.16 0.96 0.96 0.97 0.55 25 0.0033 0.0099 0.56
## BDI.17 0.96 0.96 0.97 0.55 24 0.0034 0.0106 0.55
## BDI.18 0.96 0.96 0.97 0.55 25 0.0033 0.0099 0.56
## BDI.19 0.96 0.96 0.97 0.54 24 0.0034 0.0106 0.54
## BDI.20 0.96 0.96 0.97 0.54 23 0.0035 0.0107 0.54
## BDI.21 0.96 0.96 0.97 0.55 25 0.0033 0.0096 0.56
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## BDI.1 275 0.73 0.74 0.72 0.71 0.56 0.80
## BDI.2 275 0.80 0.79 0.78 0.77 0.85 0.95
## BDI.3 275 0.84 0.83 0.83 0.81 0.93 1.03
## BDI.4 275 0.84 0.84 0.83 0.82 0.71 0.84
## BDI.5 275 0.73 0.73 0.71 0.70 0.60 0.77
## BDI.6 275 0.72 0.72 0.70 0.69 0.56 0.95
## BDI.7 275 0.84 0.83 0.83 0.81 0.88 1.07
## BDI.8 275 0.84 0.83 0.83 0.81 0.76 0.94
## BDI.9 275 0.70 0.70 0.69 0.67 0.37 0.66
## BDI.10 275 0.68 0.68 0.65 0.64 0.47 0.89
## BDI.11 275 0.73 0.74 0.73 0.71 0.56 0.74
## BDI.12 275 0.85 0.85 0.85 0.83 0.78 0.90
## BDI.13 275 0.78 0.78 0.77 0.76 0.52 0.85
## BDI.14 275 0.83 0.83 0.83 0.81 0.75 1.02
## BDI.15 275 0.80 0.80 0.79 0.77 0.89 0.92
## BDI.16 275 0.62 0.62 0.59 0.58 0.88 0.94
## BDI.17 275 0.70 0.71 0.69 0.67 0.51 0.75
## BDI.18 275 0.61 0.62 0.59 0.57 0.53 0.80
## BDI.19 275 0.74 0.74 0.73 0.71 0.65 0.84
## BDI.20 275 0.81 0.81 0.81 0.79 0.79 0.90
## BDI.21 275 0.61 0.61 0.58 0.57 0.64 0.93
##
## Non missing response frequency for each item
## 0 0.5 1 1.05 1.35 2 3 miss
## BDI.1 0.58 0 0.32 0 0 0.05 0.05 0
## BDI.2 0.45 0 0.34 0 0 0.13 0.08 0
## BDI.3 0.46 0 0.26 0 0 0.17 0.11 0
## BDI.4 0.51 0 0.30 0 0 0.16 0.03 0
## BDI.5 0.56 0 0.32 0 0 0.10 0.03 0
## BDI.6 0.68 0 0.17 0 0 0.07 0.09 0
## BDI.7 0.52 0 0.21 0 0 0.16 0.12 0
## BDI.8 0.53 0 0.24 0 0 0.17 0.06 0
## BDI.9 0.71 0 0.24 0 0 0.03 0.02 0
## BDI.10 0.72 0 0.17 0 0 0.03 0.08 0
## BDI.11 0.57 0 0.32 0 0 0.09 0.02 0
## BDI.12 0.48 0 0.30 0 0 0.16 0.05 0
## BDI.13 0.65 0 0.23 0 0 0.06 0.06 0
## BDI.14 0.58 0 0.19 0 0 0.13 0.10 0
## BDI.15 0.42 0 0.32 0 0 0.19 0.06 0
## BDI.16 0.43 0 0.32 0 0 0.17 0.07 0
## BDI.17 0.61 0 0.28 0 0 0.08 0.03 0
## BDI.18 0.63 0 0.25 0 0 0.08 0.04 0
## BDI.19 0.56 0 0.27 0 0 0.14 0.03 0
## BDI.20 0.47 0 0.33 0 0 0.13 0.06 0
## BDI.21 0.61 0 0.20 0 0 0.12 0.07 0psych::alpha(qdf[, grep("GAD", names(qdf))])##
## Reliability analysis
## Call: psych::alpha(x = qdf[, grep("GAD", names(qdf))])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.94 0.94 0.94 0.67 15 0.0058 0.84 0.83 0.66
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.92 0.94 0.95
## Duhachek 0.92 0.94 0.95
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## GAD_1 0.92 0.92 0.92 0.66 11 0.0074 0.0074 0.65
## GAD_2 0.92 0.92 0.91 0.66 11 0.0074 0.0059 0.65
## GAD_3 0.92 0.92 0.91 0.65 11 0.0075 0.0054 0.62
## GAD_4 0.93 0.93 0.92 0.67 12 0.0068 0.0110 0.62
## GAD_5 0.93 0.93 0.93 0.70 14 0.0060 0.0081 0.71
## GAD_6 0.93 0.93 0.93 0.70 14 0.0060 0.0084 0.71
## GAD_7 0.93 0.93 0.93 0.68 13 0.0067 0.0100 0.66
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## GAD_1 275 0.90 0.89 0.88 0.86 1.01 0.99
## GAD_2 275 0.90 0.89 0.89 0.86 0.87 1.02
## GAD_3 275 0.91 0.90 0.91 0.87 0.99 1.05
## GAD_4 275 0.85 0.85 0.82 0.80 0.95 1.00
## GAD_5 275 0.77 0.78 0.73 0.70 0.51 0.88
## GAD_6 275 0.78 0.78 0.72 0.70 0.83 0.93
## GAD_7 275 0.84 0.84 0.80 0.78 0.71 0.93
##
## Non missing response frequency for each item
## 0 1 1.66666666666667 2 3 miss
## GAD_1 0.37 0.36 0 0.16 0.11 0
## GAD_2 0.47 0.30 0 0.11 0.12 0
## GAD_3 0.43 0.29 0 0.16 0.13 0
## GAD_4 0.41 0.33 0 0.14 0.11 0
## GAD_5 0.69 0.16 0 0.09 0.06 0
## GAD_6 0.45 0.33 0 0.13 0.08 0
## GAD_7 0.55 0.28 0 0.11 0.07 0psych::alpha(qdf[, grep("RRS", names(qdf))])##
## Reliability analysis
## Call: psych::alpha(x = qdf[, grep("RRS", names(qdf))])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.9 0.9 0.92 0.46 8.6 0.0088 2 0.71 0.48
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.88 0.9 0.92
## Duhachek 0.88 0.9 0.92
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## RRS_1 0.89 0.89 0.91 0.47 8.1 0.0093 0.029 0.50
## RRS_2 0.88 0.88 0.90 0.45 7.2 0.0101 0.029 0.44
## RRS_3 0.88 0.88 0.90 0.45 7.4 0.0101 0.028 0.44
## RRS_4 0.89 0.88 0.90 0.45 7.5 0.0098 0.028 0.48
## RRS_5 0.91 0.91 0.92 0.52 9.7 0.0085 0.014 0.51
## RRS_6 0.89 0.89 0.91 0.46 7.7 0.0097 0.028 0.47
## RRS_7 0.88 0.88 0.90 0.45 7.4 0.0101 0.026 0.47
## RRS_8 0.89 0.88 0.90 0.45 7.5 0.0100 0.026 0.46
## RRS_9 0.88 0.88 0.90 0.44 7.2 0.0103 0.029 0.44
## RRS_10 0.89 0.89 0.90 0.46 7.8 0.0094 0.030 0.49
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## RRS_1 275 0.67 0.66 0.62 0.58 1.8 0.94
## RRS_2 275 0.79 0.80 0.78 0.73 2.2 0.98
## RRS_3 275 0.78 0.77 0.75 0.72 2.0 0.98
## RRS_4 275 0.75 0.76 0.75 0.68 2.2 0.95
## RRS_5 275 0.39 0.43 0.33 0.30 1.4 0.70
## RRS_6 275 0.72 0.71 0.67 0.64 2.4 1.02
## RRS_7 275 0.79 0.77 0.76 0.72 2.0 1.07
## RRS_8 275 0.78 0.76 0.74 0.70 2.1 1.08
## RRS_9 275 0.81 0.81 0.79 0.75 2.0 1.03
## RRS_10 275 0.69 0.71 0.68 0.61 2.2 0.96
##
## Non missing response frequency for each item
## 1 1.33333333333333 2 2.22222222222222 3 4 miss
## RRS_1 0.51 0 0.30 0 0.11 0.08 0
## RRS_2 0.27 0 0.37 0 0.24 0.12 0
## RRS_3 0.40 0 0.32 0 0.19 0.09 0
## RRS_4 0.28 0 0.40 0 0.22 0.11 0
## RRS_5 0.67 0 0.27 0 0.04 0.03 0
## RRS_6 0.21 0 0.35 0 0.25 0.19 0
## RRS_7 0.44 0 0.27 0 0.15 0.14 0
## RRS_8 0.40 0 0.27 0 0.19 0.15 0
## RRS_9 0.40 0 0.29 0 0.18 0.12 0
## RRS_10 0.28 0 0.36 0 0.25 0.11 0Sum scores
qdf$BDI_sum <- rowSums(qdf[, grep("BDI", names(qdf))])
qdf$GAD_sum <- rowSums(qdf[, grep("GAD", names(qdf))])
qdf$RRS_sum <- rowSums(qdf[, grep("RRS", names(qdf))])
qdf$depr_anx_sum <- qdf$BDI_sum + qdf$GAD_sum#any(is.na(qdf))
#write.csv(qdf, "../data/questionnaire_df_with_impute.csv")Put depr/anxisety sum and GAD sum into learn df
unique_ids <- unique(learn_df$ID)
for (i in 1:length(unique_ids)) {
learn_df[learn_df$ID==unique_ids[i], "depr_anx_sum"] <-
qdf[qdf$ID==unique_ids[i], "depr_anx_sum"]
learn_df[learn_df$ID==unique_ids[i], "bdi_sum"] <- qdf[qdf$ID==unique_ids[i], "BDI_sum"]
learn_df[learn_df$ID==unique_ids[i], "gad_sum"] <- qdf[qdf$ID==unique_ids[i], "GAD_sum"]
learn_df[learn_df$ID==unique_ids[i], "rrs_sum"] <- qdf[qdf$ID==unique_ids[i], "RRS_sum"]
test_df[test_df$ID==unique_ids[i], "depr_anx_sum"] <-
qdf[qdf$ID==unique_ids[i], "depr_anx_sum"]
test_df[test_df$ID==unique_ids[i], "bdi_sum"] <- qdf[qdf$ID==unique_ids[i], "BDI_sum"]
test_df[test_df$ID==unique_ids[i], "gad_sum"] <- qdf[qdf$ID==unique_ids[i], "GAD_sum"]
test_df[test_df$ID==unique_ids[i], "rrs_sum"] <- qdf[qdf$ID==unique_ids[i], "RRS_sum"]
}# Check
# hist(learn_df$depr_anx_sum)
# hist(test_df$depr_anx_sum)
#
# hist(learn_df$rrs_sum)
# hist(test_df$rrs_sum)# Export — df imported into for brms_s
# write.csv(learn_df, "../data/learn_df_with_qairre_data.csv")
# write.csv(test_df, "../data/test_df_with_qairre_data.csv")Plot for supp fig
bdi_means <- qdf %>% group_by(group) %>% summarize(m=mean(BDI_sum))
bdi_p <- ggplot(qdf, aes(x=BDI_sum, fill=group, color=group)) +
geom_vline(data=bdi_means, aes(xintercept = m, color=group), size=3) +
geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
xlab("") +
theme(axis.text.y = element_blank()) +
theme(axis.ticks.y = element_blank()) +
scale_fill_manual(values=c("red", "blue", "gray")) +
scale_color_manual(values=c("red", "blue", "gray")) +
ggtitle("Depression (BDI-II)") + tp +
tolJust for legend
# ggplot(qdf, aes(x=BDI_sum, fill=group, color=group)) +
# geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
# xlab("") +
# theme(axis.text.y = element_blank()) +
# theme(axis.ticks.y = element_blank()) +
# scale_fill_manual(values=c("red", "blue", "gray")) +
# scale_color_manual(values=c("red", "blue", "gray")) +
# ggtitle("Depression (BDI-II)") + tpgad_means <- qdf %>% group_by(group) %>% summarize(m=mean(GAD_sum))
gad_p <- ggplot(qdf, aes(x=GAD_sum, fill=group, color=group)) +
geom_vline(data=gad_means, aes(xintercept = m, color=group), size=3) +
geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
xlab("") +
theme(axis.text.y = element_blank()) +
theme(axis.ticks.y = element_blank()) +
scale_fill_manual(values=c("red", "blue", "gray")) +
scale_color_manual(values=c("red", "blue", "gray")) +
ggtitle("Generalized anxiety (GAD-7)") + tp +
tolrrs_means <- qdf %>% group_by(group) %>% summarize(m=mean(RRS_sum))
rrs_p <- ggplot(qdf, aes(x=RRS_sum, fill=group, color=group)) +
geom_vline(data=rrs_means, aes(xintercept = m, color=group), size=3) +
geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
xlab("") +
theme(axis.text.y = element_blank()) +
theme(axis.ticks.y = element_blank()) +
scale_fill_manual(values=c("red", "blue", "gray")) +
scale_color_manual(values=c("red", "blue", "gray")) +
ggtitle("Rumination (RRS-SF)") + tp +
tolsum(table(qdf$group))## [1] 275sym_plot <- bdi_p/ gad_p / rrs_psym_plot
#ggsave("../paper/figs/supp-figs/symptom_plot.png", sym_plot, width=8.5, height = 8.5, dpi=200)Behavioral performance and individual differences
cor.test(learn_m$m, qdf$depr_anx_sum)##
## Pearson's product-moment correlation
##
## data: learn_m$m and qdf$depr_anx_sum
## t = -0.01433, df = 273, p-value = 0.9886
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1191390 0.1174287
## sample estimates:
## cor
## -0.0008672926cor.test(test_m$m, qdf$depr_anx_sum)##
## Pearson's product-moment correlation
##
## data: test_m$m and qdf$depr_anx_sum
## t = 0.14914, df = 273, p-value = 0.8816
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1093747 0.1271742
## sample estimates:
## cor
## 0.009026022cor.test(learn_m$m, qdf$RRS_sum)##
## Pearson's product-moment correlation
##
## data: learn_m$m and qdf$RRS_sum
## t = -0.023301, df = 273, p-value = 0.9814
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1196743 0.1168932
## sample estimates:
## cor
## -0.001410259cor.test(test_m$m, qdf$RRS_sum)##
## Pearson's product-moment correlation
##
## data: test_m$m and qdf$RRS_sum
## t = 0.0072235, df = 273, p-value = 0.9942
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1178529 0.1187150
## sample estimates:
## cor
## 0.0004371846ctest <- cor.test(qdf$GAD_sum, qdf$BDI_sum)
r <- round(ctest$estimate, 2)
p <- round(ctest$p.value, 2)
str <- paste("r =", r, "p =", p)
ggplot(qdf, aes(x=GAD_sum, BDI_sum)) +
geom_smooth(method='lm', formula= y~x, color="black", size=3) +
geom_point(size=6, fill="white", pch=21) +
ga + ap + tp +
ggtitle("", subtitle=str) +
ylab("BDI-II") + xlab("GAD-7") +
theme(plot.subtitle = element_text(size = 25, face = "bold"))
Create high low DA, BDI and GAD
median(qdf$BDI_sum)## [1] 10median(qdf$GAD_sum)## [1] 4hist(qdf$depr_anx_sum, breaks=100)
median(qdf$depr_anx_sum)## [1] 15table(if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0))##
## 0 1
## 138 137#table(if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0)) # Check that's about half
qdf$DA_hilo <- if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0)
qdf %>% group_by(DA_hilo) %>% summarize(m=median(depr_anx_sum)) # Confirm## # A tibble: 2 × 2
## DA_hilo m
## <dbl> <dbl>
## 1 0 4
## 2 1 33#table(if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0))
qdf$dep_hilo <- if_else(qdf$BDI_sum > median(qdf$BDI_sum), 1, 0)
#qdf %>% group_by(dep_hilo) %>% summarize(m=median(depr_anx_sum)) #table(if_else(qdf$GAD_sum > median(qdf$GAD_sum), 1, 0))
qdf$anx_hilo <- if_else(qdf$GAD_sum > median(qdf$GAD_sum), 1, 0)
#qdf %>% group_by(anx_hilo) %>% summarize(m=median(GAD_sum)) table(if_else(qdf$RRS_sum > median(qdf$RRS_sum), 1, 0))##
## 0 1
## 146 129qdf$rrs_hilo <- if_else(qdf$RRS_sum > median(qdf$RRS_sum), 1, 0)
qdf %>% group_by(rrs_hilo) %>% summarize(m=median(RRS_sum)) ## # A tibble: 2 × 2
## rrs_hilo m
## <dbl> <dbl>
## 1 0 15
## 2 1 25Put into learning and test dfs
high_depr_anx_pts <- qdf[qdf$DA_hilo == 1, "ID"]
low_depr_anx_pts <- qdf[qdf$DA_hilo == 0, "ID"]
learn_df[learn_df$ID %in% high_depr_anx_pts, "DA_hilo"] <- "high depr-anx"
learn_df[learn_df$ID %in% low_depr_anx_pts, "DA_hilo"] <- "low depr-anx"
test_df[test_df$ID %in% high_depr_anx_pts, "DA_hilo"] <- "high depr-anx"
test_df[test_df$ID %in% low_depr_anx_pts, "DA_hilo"] <- "low depr-anx"high_rum_pts <- qdf[qdf$rrs_hilo == 1, "ID"]
low_rum_pts <- qdf[qdf$rrs_hilo == 0, "ID"]
learn_df[learn_df$ID %in% high_rum_pts, "rrs_hilo"] <- "high rum"
learn_df[learn_df$ID %in% low_rum_pts, "rrs_hilo"] <- "low rum"
test_df[test_df$ID %in% high_rum_pts, "rrs_hilo"] <- "high rum"
test_df[test_df$ID %in% low_rum_pts, "rrs_hilo"] <- "low rum"Plot proportion correct by median split
pcor_ss_DA <- data.frame(learn_df %>% group_by(stim_iter, set_size,
ID, DA_hilo) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'stim_iter', 'set_size', 'ID'. You can
## override using the `.groups` argument.pcor_ss_err_DA <- Rmisc::summarySEwithin(pcor_ss_DA,
measurevar = "m",
withinvars = c("set_size", "stim_iter", "DA_hilo"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_size, stim_iter, DA_hilopcor_ss_err_DA$DA_hilo <- factor(pcor_ss_err_DA$DA_hilo, levels=c("low depr-anx", "high depr-anx"))pcor_by_grp <- ggplot(pcor_ss_err_DA, aes(x=stim_iter, y=m, group=as.factor(DA_hilo),
color=as.factor(DA_hilo))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_point(aes(fill=as.factor(DA_hilo)), color="black", size=5, pch=21) +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
#geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") +
ggtitle("") + tp + facet_wrap(~ set_size) + ft + lp +
scale_fill_manual(values=c("orange1", "brown")) +
scale_color_manual(values=c("orange1", "brown")) + tolpcor_by_grp 
#ggsave("../paper/figs/pieces/pcor_by_grp.png", pcor_by_grp, width=14, height=10, dpi=300)pcor_ss_rrs <- data.frame(learn_df %>% group_by(stim_iter, set_size,
ID, rrs_hilo) %>% summarize(m=mean(correct)))## `summarise()` has grouped output by 'stim_iter', 'set_size', 'ID'. You can
## override using the `.groups` argument.pcor_ss_err_rrs <- Rmisc::summarySEwithin(pcor_ss_rrs,
measurevar = "m",
withinvars = c("set_size", "stim_iter", "rrs_hilo"),
idvar = "ID")## Automatically converting the following non-factors to factors: set_size, stim_iter, rrs_hilopcor_ss_err_rrs$rrs_hilo <- factor(pcor_ss_err_rrs$rrs_hilo, levels=c("low rum", "high rum"))pcor_by_grp_rum <- ggplot(pcor_ss_err_rrs, aes(x=stim_iter, y=m, group=as.factor(rrs_hilo),
color=as.factor(rrs_hilo))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_point(aes(fill=as.factor(rrs_hilo)), color="black", size=5, pch=21) +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
#geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") +
ggtitle("") + tp + facet_wrap(~ set_size) + ft + lp +
scale_fill_manual(values=c("red", "darkred")) +
scale_color_manual(values=c("red", "darkred"))Rumination version
pcor_by_grp_rum
#ggsave("../paper/figs/pieces/pcor_by_grp_rum.png", pcor_by_grp_rum, width=14, height=8, dpi=300)test_error_df <- test_df %>% filter(correct==0)
np_test_summs_da <- test_error_df %>% group_by(ID, DA_hilo, phase) %>%
summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))## `summarise()` has grouped output by 'ID', 'DA_hilo'. You can override using the
## `.groups` argument.np_test_si_da <- Rmisc::summarySEwithin(np_test_summs_da,
measurevar = "neutral_pref",
withinvars = c("DA_hilo", "phase"),
idvar = "ID")## Automatically converting the following non-factors to factors: DA_hilo, phasenp_test_si_da$DA_hilo <- factor(np_test_si_da$DA_hilo, levels=c("low depr-anx", "high depr-anx"))ggplot(np_test_si_da, aes(x=phase, y=neutral_pref, fill=phase)) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=neutral_pref-se, ymax=neutral_pref+se), width=.2) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(#"Neutral preference during test
"Empirical") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23) + facet_wrap(~ DA_hilo) + ft
Read in BRMS posteriors fit with brms_s.R for
behavior
Performace in terms of proprortion correct
phase1_learn_perf <- read.csv("../model_res/brms_res/reduced_perf_model_phase1__86050__.csv")
phase2_learn_perf <- read.csv("../model_res/brms_res/reduced_correct_perf_model_phase2__76459__.csv")
test_perf <- read.csv("../model_res/brms_res/reduced_correct_test_model__24237__.csv")ReturnPosteriorMeanAnd90CI <- function(posterior) {
cat("\nposterior mean =", mean(posterior))
cat("\n90% CI ="); print (bayestestR::ci(posterior, ci = .9, method = "HDI"))
}Neutral preference
np_da <- read.csv("../model_res/brms_res/neutral_pref_learn_model__64424__.csv")
# Version that converged that has intercept set to 0 but models REs in depr/anx
np_test_da <- read.csv("../model_res/brms_res/neutral_pref_test_model_0int__14655__.csv")Rumination
phase1_rum <- read.csv("../model_res/brms_res/RUM_reduced_perf_model_phase1__63279__.csv")
phase2_rum <- read.csv("../model_res/brms_res/RUM_reduced_correct_perf_model_phase2__25543__.csv")
test_rum <- read.csv("../model_res/brms_res/RUM_reduced_correct_test_model__53350__.csv")np_rum <- read.csv("../model_res/brms_res/RUM_neutral_pref_learn_model__89767__.csv")
np_test_rum <- read.csv("../model_res/brms_res/RUM_neutral_pref_test_model__71650__.csv")Proportion correct ~ symptoms
Phase 1
cat("\n Depression prop correct phase 1: \n")##
## Depression prop correct phase 1:# DA
ReturnPosteriorMeanAnd90CI(phase1_learn_perf$b_scaledepr_anx_sum)##
## posterior mean = -0.005247766
## 90% CI =90% HDI: [-0.04, 0.03]ReturnPosteriorMeanAnd90CI(phase1_learn_perf$b_scaleset_size) # Sanity check this is negative and tight ##
## posterior mean = -0.3166876
## 90% CI =90% HDI: [-0.33, -0.30]ReturnPosteriorMeanAnd90CI(phase1_learn_perf$b_scaledepr_anx_sum.scaleset_size)##
## posterior mean = 0.004837842
## 90% CI =90% HDI: [-0.01, 0.02]cat("\n\n Rumination prop correct phase 1: \n")##
##
## Rumination prop correct phase 1:# Rum
ReturnPosteriorMeanAnd90CI(phase1_rum$b_scalerrs_sum)##
## posterior mean = -0.007088953
## 90% CI =90% HDI: [-0.04, 0.03]ReturnPosteriorMeanAnd90CI(phase1_rum$b_scaleset_size) # Sanity check - v consistent with above ##
## posterior mean = -0.3169965
## 90% CI =90% HDI: [-0.33, -0.30]ReturnPosteriorMeanAnd90CI(phase1_rum$b_scalerrs_sum.scaleset_size)##
## posterior mean = 0.01584776
## 90% CI =90% HDI: [ 0.00, 0.03]Phase 2
# DA
cat("\n\n Depression prop correct phase 2: \n")##
##
## Depression prop correct phase 2:ReturnPosteriorMeanAnd90CI(phase2_learn_perf$b_scaledepr_anx_sum)##
## posterior mean = 6.641381e-05
## 90% CI =90% HDI: [-0.08, 0.07]ReturnPosteriorMeanAnd90CI(phase2_learn_perf$b_scaleset_size) # Sanity check this is negative and tight ##
## posterior mean = -0.2150235
## 90% CI =90% HDI: [-0.24, -0.19]ReturnPosteriorMeanAnd90CI(phase2_learn_perf$b_scaledepr_anx_sum.scaleset_size)##
## posterior mean = 0.02980973
## 90% CI =90% HDI: [0.01, 0.05]# Rum
cat("\n\n Rumination prop correct phase 2: \n")##
##
## Rumination prop correct phase 2:ReturnPosteriorMeanAnd90CI(phase2_rum$b_scalerrs_sum)##
## posterior mean = 0.0003879612
## 90% CI =90% HDI: [-0.08, 0.08]ReturnPosteriorMeanAnd90CI(phase2_rum$b_scaleset_size) # Sanity check - v consistent with above ##
## posterior mean = -0.2149793
## 90% CI =90% HDI: [-0.24, -0.19]ReturnPosteriorMeanAnd90CI(phase2_rum$b_scalerrs_sum.scaleset_size)##
## posterior mean = 0.01623668
## 90% CI =90% HDI: [-0.01, 0.04]Test
# DA
cat("\n\n Depression prop correct test: \n")##
##
## Depression prop correct test:ReturnPosteriorMeanAnd90CI(test_perf$b_scaledepr_anx_sum)##
## posterior mean = 0.009176028
## 90% CI =90% HDI: [-0.08, 0.09]# Rum
cat("\n\n Rumination prop correct test: \n")##
##
## Rumination prop correct test:ReturnPosteriorMeanAnd90CI(test_rum$b_scalerrs_sum)##
## posterior mean = 0.004301401
## 90% CI =90% HDI: [-0.08, 0.09]Neutral preference ~ symptoms
Learning
The intercept posteriors show that the evidence for a neutral preference found in the frequentist models holds in these Bayesian regression models, now while also controlling for psychiatric symptoms
cat("\n Depression : \n")##
## Depression :cat("\n Evidence controlling for depression : "); ReturnPosteriorMeanAnd90CI(np_da$b_Intercept)##
## Evidence controlling for depression :##
## posterior mean = 0.08774017
## 90% CI =90% HDI: [0.06, 0.11]cat("\n Evidence for depession effect : ");ReturnPosteriorMeanAnd90CI(np_da$b_scaledepr_anx_sum)##
## Evidence for depession effect :##
## posterior mean = -0.01066593
## 90% CI =90% HDI: [-0.03, 0.01]cat("\n Rumination : \n")##
## Rumination :cat("\n Evidence controlling for rumination : "); ReturnPosteriorMeanAnd90CI(np_rum$b_Intercept)##
## Evidence controlling for rumination :##
## posterior mean = 0.08747077
## 90% CI =90% HDI: [0.06, 0.11]cat("\n Evidence for rumination effect : ");ReturnPosteriorMeanAnd90CI(np_rum$b_scalerrs_sum)##
## Evidence for rumination effect :##
## posterior mean = 0.004026421
## 90% CI =90% HDI: [-0.02, 0.03]Test
Note that HDI for intercept (indicating neutral pref at test) includes 0 (although it’s close controlling for depression symptoms)
cat("\n Depression : \n")##
## Depression :cat("\n Evidence controlling for depression : "); ReturnPosteriorMeanAnd90CI(np_test_da$b_Intercept)##
## Evidence controlling for depression :##
## posterior mean = 0.03915869
## 90% CI =90% HDI: [ 0.00, 0.09]cat("\n Evidence for depession effect : ");ReturnPosteriorMeanAnd90CI(np_test_da$b_scaledepr_anx_sum)##
## Evidence for depession effect :##
## posterior mean = -0.007424603
## 90% CI =90% HDI: [-0.06, 0.05]cat("\n Rumination : \n")##
## Rumination :cat("\n Evidence controlling for rumination : "); ReturnPosteriorMeanAnd90CI(np_test_rum$b_Intercept)##
## Evidence controlling for rumination :##
## posterior mean = 0.03967129
## 90% CI =90% HDI: [-0.01, 0.09]cat("\n Evidence for rumination effect : ");ReturnPosteriorMeanAnd90CI(np_test_rum$b_scalerrs_sum)##
## Evidence for rumination effect :##
## posterior mean = -0.00610522
## 90% CI =90% HDI: [-0.06, 0.05]Make sure Rhats below 1.1
CheckRhat("../model_res/brms_res/Rhatreduced_correct_perf_model_phase2__76459__.csv")## [1] TRUECheckRhat("../model_res/brms_res/Rhatreduced_perf_model_phase1__86050__.csv")## [1] TRUECheckRhat("../model_res/brms_res/Rhatreduced_correct_test_model__24237__.csv")## [1] TRUECheckRhat("../model_res/brms_res/Rhatneutral_pref_learn_model__64424__.csv")## [1] TRUECheckRhat("../model_res/brms_res/Rhatneutral_pref_test_model__70753__.csv")## [1] TRUECheckRhat("../model_res/brms_res/RhatRUM_neutral_pref_learn_model__89767__.csv")## [1] TRUECheckRhat("../model_res/brms_res/RhatRUM_neutral_pref_test_model__71650__.csv")## [1] TRUECheckRhat("../model_res/brms_res/RhatRUM_reduced_correct_perf_model_phase2__25543__.csv")## [1] TRUECheckRhat("../model_res/brms_res/RhatRUM_reduced_perf_model_phase1__63279__.csv")## [1] TRUECheckRhat("../model_res/brms_res/RhatRUM_reduced_correct_test_model__53350__.csv")## [1] TRUEGet ready to plot proportion correct posterior for depr/anxiety …
depr_anx_learn_perf_posteriors <- rbind(
data.frame("trace"=phase1_learn_perf$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_learn_perf$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Learn Phase 2"),
data.frame("trace"=phase1_learn_perf$b_scaledepr_anx_sum.scaleset_size,
"label"="da_ss_int", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_learn_perf$b_scaledepr_anx_sum.scaleset_size,
"label"="da_ss_int", "phase"="Learn Phase 2")
)
depr_anx_test_perf_posterior <- rbind(
data.frame("trace"=test_perf$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Test")
) learn_correct_posterior_p <- ggplot(depr_anx_learn_perf_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("tan1", "chocolate")) + ylab("") +
theme(axis.text.x = element_text(size=15), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
learn_correct_posterior_p
test_correct_posterior_p <-ggplot(depr_anx_test_perf_posterior,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) +
ga + ap +
scale_color_manual(values=c("tan1")) +
ylab("") +
theme(axis.text.x = element_text(size=15),
axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Test") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp #+
test_correct_posterior_p 
# ggsave("../paper/figs/pieces/fig7_learn_correct_posterior_p.png", learn_correct_posterior_p, width=2.2, height=3.3, dpi=500)
# ggsave("../paper/figs/pieces/fig7_test_correct_posterior_p.png",
# test_correct_posterior_p, width=2.2, height=1.65, dpi=500)… and rum
rum_learn_perf_posteriors <- rbind(
data.frame("trace"=phase1_rum$b_scalerrs_sum, "label"="rum_main_effect", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_rum$b_scalerrs_sum, "label"="rum_main_effect", "phase"="Learn Phase 2"),
data.frame("trace"=phase1_rum$b_scalerrs_sum.scaleset_size,
"label"="rum_ss_int", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_rum$b_scalerrs_sum.scaleset_size,
"label"="rum_ss_int", "phase"="Learn Phase 2")
)
rum_test_perf_posterior <- rbind(
data.frame("trace"=test_rum$b_scalerrs_sum, "label"="rum_main_effect", "phase"="Test")
) learn_correct_posterior_p_RUM <-
ggplot(rum_learn_perf_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("indianred2", "darkred")) + ylab("") +
theme(axis.text.x = element_text(size=15), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
learn_correct_posterior_p_RUM
test_correct_posterior_p_RUM <-
ggplot(rum_test_perf_posterior,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) +
ga + ap +
scale_color_manual(values=c("indianred1")) +
ylab("") +
theme(axis.text.x = element_text(size=15), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Test") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp #+
test_correct_posterior_p_RUM
# ggsave("../paper/figs/pieces/fig7_learn_correct_posterior_p_RUM.png", learn_correct_posterior_p_RUM, width=2.2, height=3.3, dpi=500)
# ggsave("../paper/figs/pieces/fig7_test_correct_posterior_p_RUM.png",
# test_correct_posterior_p_RUM, width=2.2, height=1.65, dpi=500)Get ready to plot neutral pref posterior for depr/anxiety and rum
depr_anx_np_posteriors <- rbind(
data.frame("trace"=np_da$b_scaledepr_anx_sum, "label"="da_effect", "phase"="Learn"),
data.frame("trace"=np_test_da$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Test")
)
rum_np_posteriors <- rbind(
data.frame("trace"=np_rum$b_scalerrs_sum, "label"="rum_effect", "phase"="Learn"),
data.frame("trace"=np_test_rum$b_scalerrs_sum, "label"="da_main_effect", "phase"="Test")
) np_da_posterior_p <-
ggplot(depr_anx_np_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("tan1", "tan1")) + ylab("") +
theme(axis.text.x = element_text(size=10), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
np_da_posterior_p 
rum_np_posterior_p <-
ggplot(rum_np_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("indianred2", "indianred2")) + ylab("") +
theme(axis.text.x = element_text(size=10), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
rum_np_posterior_p
# ggsave("../paper/figs/pieces/fig7_np_da_posterior_p.png", np_da_posterior_p ,
# width=2.2, height=3.3, dpi=500)
# ggsave("../paper/figs/pieces/fig7_rum_np_posterior_p.png",
# rum_np_posterior_p, width=2.2, height=3.3, dpi=500)Symptoms ~ model parameters
learn_m <- learn_df %>% group_by(ID) %>% summarize(m=mean(correct))
assert(all(learn_m$ID==qdf$ID))
test_m <- test_df %>% group_by(ID) %>% summarize(m=mean(correct))
assert(all(learn_m$ID==qdf$ID))
assert(all(m35$ID==qdf$ID))
assert(all(m35$ID==learn_m$ID))
assert(all(m35$ID==test_m$ID))
m35$depr_anx_sum <- qdf$depr_anx_sum
m35$BDI_sum <- qdf$BDI_sum
m35$GAD_sum <- qdf$GAD_sum
m35$RRS_sum <- qdf$RRS_sum
m35$learn_m <- learn_m$m
m35$test_m <- test_m$mMany parameters predict learning and test performance…
summary(learn_preds_bayes <- rstanarm::stan_glm(
scale(learn_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.034 seconds (Warm-up)
## Chain 1: 0.052 seconds (Sampling)
## Chain 1: 0.086 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 2: 0.054 seconds (Sampling)
## Chain 2: 0.086 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 9e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.031 seconds (Warm-up)
## Chain 3: 0.052 seconds (Sampling)
## Chain 3: 0.083 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 9e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 4: 0.047 seconds (Sampling)
## Chain 4: 0.08 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(learn_m) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.0 0.0 0.0 0.0
## scale(kappa) 0.3 0.0 0.2 0.3 0.3
## scale(alpha_pos) 0.4 0.0 0.4 0.4 0.5
## scale(alpha_neg) 0.2 0.0 0.1 0.2 0.2
## scale(phi) -0.3 0.0 -0.3 -0.3 -0.2
## scale(rho) 0.1 0.0 0.1 0.1 0.2
## scale(rl_off) -0.3 0.0 -0.3 -0.3 -0.2
## scale(epsilon) -0.3 0.0 -0.4 -0.3 -0.3
## scale(not_pun_bonus) -0.1 0.0 -0.2 -0.1 -0.1
## sigma 0.5 0.0 0.5 0.5 0.6
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.0 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4451
## scale(kappa) 0.0 1.0 3653
## scale(alpha_pos) 0.0 1.0 3728
## scale(alpha_neg) 0.0 1.0 4008
## scale(phi) 0.0 1.0 3418
## scale(rho) 0.0 1.0 3689
## scale(rl_off) 0.0 1.0 4010
## scale(epsilon) 0.0 1.0 3723
## scale(not_pun_bonus) 0.0 1.0 3857
## sigma 0.0 1.0 3701
## mean_PPD 0.0 1.0 4174
## log-posterior 0.1 1.0 1740
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(test_preds_bayes <- rstanarm::stan_glm(
scale(test_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 1: 0.048 seconds (Sampling)
## Chain 1: 0.08 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 2: 0.047 seconds (Sampling)
## Chain 2: 0.08 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 3: 0.044 seconds (Sampling)
## Chain 3: 0.077 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.034 seconds (Warm-up)
## Chain 4: 0.053 seconds (Sampling)
## Chain 4: 0.087 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(test_m) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.0 0.0 0.0 0.0
## scale(kappa) 0.1 0.0 0.0 0.1 0.1
## scale(alpha_pos) 0.6 0.0 0.6 0.6 0.7
## scale(alpha_neg) 0.2 0.0 0.1 0.2 0.2
## scale(phi) -0.2 0.0 -0.2 -0.2 -0.1
## scale(rho) 0.1 0.0 0.0 0.1 0.1
## scale(rl_off) -0.5 0.0 -0.5 -0.5 -0.4
## scale(epsilon) -0.2 0.0 -0.2 -0.2 -0.1
## scale(not_pun_bonus) 0.0 0.0 -0.1 0.0 0.0
## sigma 0.5 0.0 0.5 0.5 0.6
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.0 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3958
## scale(kappa) 0.0 1.0 3307
## scale(alpha_pos) 0.0 1.0 3426
## scale(alpha_neg) 0.0 1.0 4644
## scale(phi) 0.0 1.0 3455
## scale(rho) 0.0 1.0 3955
## scale(rl_off) 0.0 1.0 3959
## scale(epsilon) 0.0 1.0 4398
## scale(not_pun_bonus) 0.0 1.0 4318
## sigma 0.0 1.0 4648
## mean_PPD 0.0 1.0 4140
## log-posterior 0.1 1.0 1777
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.03 seconds (Warm-up)
## Chain 1: 0.047 seconds (Sampling)
## Chain 1: 0.077 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.031 seconds (Warm-up)
## Chain 2: 0.049 seconds (Sampling)
## Chain 2: 0.08 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.059 seconds (Warm-up)
## Chain 3: 0.05 seconds (Sampling)
## Chain 3: 0.109 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 4: 0.056 seconds (Sampling)
## Chain 4: 0.088 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(kappa) -0.1 0.1 -0.2 -0.1 0.0
## scale(alpha_pos) 0.1 0.1 0.0 0.1 0.1
## scale(alpha_neg) 0.0 0.1 -0.1 0.0 0.1
## scale(phi) -0.1 0.1 -0.2 -0.1 0.0
## scale(rho) 0.0 0.1 -0.1 0.0 0.1
## scale(rl_off) 0.0 0.1 -0.1 0.0 0.1
## scale(epsilon) 0.0 0.1 0.0 0.0 0.1
## scale(not_pun_bonus) 0.1 0.1 0.0 0.1 0.1
## sigma 1.0 0.0 1.0 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3809
## scale(kappa) 0.0 1.0 3137
## scale(alpha_pos) 0.0 1.0 3216
## scale(alpha_neg) 0.0 1.0 3906
## scale(phi) 0.0 1.0 3193
## scale(rho) 0.0 1.0 3812
## scale(rl_off) 0.0 1.0 4305
## scale(epsilon) 0.0 1.0 4081
## scale(not_pun_bonus) 0.0 1.0 4331
## sigma 0.0 1.0 3725
## mean_PPD 0.0 1.0 3910
## log-posterior 0.1 1.0 1725
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(rrs_preds_bayes <- rstanarm::stan_glm(
scale(RRS_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.031 seconds (Warm-up)
## Chain 1: 0.05 seconds (Sampling)
## Chain 1: 0.081 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.03 seconds (Warm-up)
## Chain 2: 0.047 seconds (Sampling)
## Chain 2: 0.077 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 3: 0.046 seconds (Sampling)
## Chain 3: 0.078 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 5e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.034 seconds (Warm-up)
## Chain 4: 0.05 seconds (Sampling)
## Chain 4: 0.084 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(RRS_sum) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(kappa) 0.0 0.1 -0.1 0.0 0.1
## scale(alpha_pos) 0.0 0.1 -0.1 0.0 0.1
## scale(alpha_neg) 0.0 0.1 -0.1 0.0 0.1
## scale(phi) 0.0 0.1 -0.1 0.0 0.1
## scale(rho) 0.0 0.1 0.0 0.0 0.1
## scale(rl_off) 0.0 0.1 0.0 0.0 0.1
## scale(epsilon) 0.0 0.1 0.0 0.0 0.1
## scale(not_pun_bonus) 0.0 0.1 -0.1 0.0 0.1
## sigma 1.0 0.0 1.0 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4584
## scale(kappa) 0.0 1.0 3550
## scale(alpha_pos) 0.0 1.0 2929
## scale(alpha_neg) 0.0 1.0 4164
## scale(phi) 0.0 1.0 3188
## scale(rho) 0.0 1.0 3853
## scale(rl_off) 0.0 1.0 3976
## scale(epsilon) 0.0 1.0 4285
## scale(not_pun_bonus) 0.0 1.0 3743
## sigma 0.0 1.0 3779
## mean_PPD 0.0 1.0 4241
## log-posterior 0.1 1.0 1875
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).packageVersion("rstanarm")## [1] '2.32.1'Get traces for each par
learn_traces <- data.frame(as.matrix(learn_preds_bayes))
learn_est <- ExtractEstimates(learn_traces)
learn_est$label <- "Learning"
learn_est$category <- "Proportion correct"
# Drop the coefs that aren't model pars
learn_est <- learn_est %>% filter(!coef %in% c("sigma", "X.Intercept."))test_traces <- data.frame(as.matrix(test_preds_bayes))
test_est <- ExtractEstimates(test_traces)
test_est$label <- "Test"
test_est$category <- "Proportion correct"
test_est <- test_est %>% filter(!coef %in% c("sigma", "X.Intercept."))
assert(all(test_est$coef==learn_est$coef))da_traces <- data.frame(as.matrix(da_preds_bayes))
da_est <- ExtractEstimates(da_traces)
da_est$label <- "Depr/anx"
da_est$category <- "Symptoms"
da_est <- da_est %>% filter(!coef %in% c("sigma", "X.Intercept."))
assert(all(da_est$coef==learn_est$coef))ReturnPosteriorMeanAnd90CI(da_traces$scale.alpha_neg.)##
## posterior mean = 0.005131044
## 90% CI =90% HDI: [-0.09, 0.11]ReturnPosteriorMeanAnd90CI(da_traces$scale.alpha_pos.) ##
## posterior mean = 0.05431283
## 90% CI =90% HDI: [-0.06, 0.17]ReturnPosteriorMeanAnd90CI(da_traces$scale.kappa.)##
## posterior mean = -0.08965521
## 90% CI =90% HDI: [-0.21, 0.03]ReturnPosteriorMeanAnd90CI(da_traces$scale.phi.)##
## posterior mean = -0.07042382
## 90% CI =90% HDI: [-0.18, 0.05]ReturnPosteriorMeanAnd90CI(da_traces$scale.rho.)##
## posterior mean = 0.008827428
## 90% CI =90% HDI: [-0.10, 0.11]ReturnPosteriorMeanAnd90CI(da_traces$scale.rl_off.)##
## posterior mean = 0.005668033
## 90% CI =90% HDI: [-0.10, 0.12]ReturnPosteriorMeanAnd90CI(da_traces$scale.epsilon.)##
## posterior mean = 0.04158132
## 90% CI =90% HDI: [-0.06, 0.15]ReturnPosteriorMeanAnd90CI(da_traces$scale.not_pun_bonus.)##
## posterior mean = 0.06771418
## 90% CI =90% HDI: [-0.04, 0.17]The three with close to some support
cat("\n Depression and bonus for non-punishing items\n")##
## Depression and bonus for non-punishing itemsReturnPosteriorMeanAnd90CI(da_traces$scale.not_pun_bonus.)##
## posterior mean = 0.06771418
## 90% CI =90% HDI: [-0.04, 0.17]cat("\n Depression and WM decay\n")##
## Depression and WM decayReturnPosteriorMeanAnd90CI(da_traces$scale.phi.)##
## posterior mean = -0.07042382
## 90% CI =90% HDI: [-0.18, 0.05]cat("\n Depression and WM capacity\n")##
## Depression and WM capacityReturnPosteriorMeanAnd90CI(da_traces$scale.kappa.)##
## posterior mean = -0.08965521
## 90% CI =90% HDI: [-0.21, 0.03]cat("\n Depression and learning rate from negative PEs \n")##
## Depression and learning rate from negative PEsReturnPosteriorMeanAnd90CI(da_traces$scale.alpha_neg.)##
## posterior mean = 0.005131044
## 90% CI =90% HDI: [-0.09, 0.11]Robustness check in univariate models the results aligned w hypotheses
summary(npb_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(not_pun_bonus),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 1: 0.034 seconds (Sampling)
## Chain 1: 0.051 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 2: 0.033 seconds (Sampling)
## Chain 2: 0.049 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 3: 0.033 seconds (Sampling)
## Chain 3: 0.049 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 4: 0.032 seconds (Sampling)
## Chain 4: 0.048 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(not_pun_bonus) 0.0 0.1 0.0 0.0 0.1
## sigma 1.0 0.0 1.0 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3933
## scale(not_pun_bonus) 0.0 1.0 4136
## sigma 0.0 1.0 3796
## mean_PPD 0.0 1.0 3975
## log-posterior 0.0 1.0 1797
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).Of note, results here and elsewhere very slightly discrepant when run within notebook vs. knitted results presumably due to rounding during knit/posterior sampling randomness but does not change results
cat("\n Univariate depression/anxiety bonus for non-punishing items\n")##
## Univariate depression/anxiety bonus for non-punishing itemsnpb_traces <- data.frame(as.matrix(npb_da_preds_bayes))
ReturnPosteriorMeanAnd90CI(npb_traces$scale.not_pun_bonus.)##
## posterior mean = 0.04539584
## 90% CI =90% HDI: [-0.06, 0.15]rrs_traces <- data.frame(as.matrix(rrs_preds_bayes))
rrs_est <- ExtractEstimates(rrs_traces)
rrs_est$label <- "Rumination"
rrs_est$category <- "Symptoms"
rrs_est <- rrs_est %>% filter(!coef %in% c("sigma", "X.Intercept."))
assert(all(rrs_est$coef==test_est$coef))ReturnPosteriorMeanAnd90CI(rrs_traces$scale.alpha_neg.)##
## posterior mean = -0.002717084
## 90% CI =90% HDI: [-0.10, 0.11]ReturnPosteriorMeanAnd90CI(rrs_traces$scale.alpha_pos.) ##
## posterior mean = 0.01449482
## 90% CI =90% HDI: [-0.09, 0.14]ReturnPosteriorMeanAnd90CI(rrs_traces$scale.kappa.)##
## posterior mean = -0.003328104
## 90% CI =90% HDI: [-0.13, 0.11]ReturnPosteriorMeanAnd90CI(rrs_traces$scale.phi.)##
## posterior mean = -0.02781697
## 90% CI =90% HDI: [-0.15, 0.09]ReturnPosteriorMeanAnd90CI(rrs_traces$scale.rho.)##
## posterior mean = 0.03524774
## 90% CI =90% HDI: [-0.07, 0.14]ReturnPosteriorMeanAnd90CI(rrs_traces$scale.rl_off.)##
## posterior mean = 0.04702406
## 90% CI =90% HDI: [-0.06, 0.16]ReturnPosteriorMeanAnd90CI(rrs_traces$scale.epsilon.)##
## posterior mean = 0.03527464
## 90% CI =90% HDI: [-0.07, 0.14]ReturnPosteriorMeanAnd90CI(rrs_traces$scale.not_pun_bonus.)##
## posterior mean = -0.001529133
## 90% CI =90% HDI: [-0.11, 0.10]all_ests <- rbind(learn_est, test_est, da_est, rrs_est)
all_ests$label <- factor(all_ests$label, levels=c("Learning", "Test", "Depr/anx", "Rumination"))all_ests_p <- ggplot(all_ests, #%>% filter(category=="Proportion correct"),
aes(y=coef, x=m, fill=label)) +
geom_vline(xintercept=0, size=1.3, color="gray35") +
geom_vline(xintercept=c(-.5, -.25, .25, .5), color="gray65") +
geom_errorbar(aes(xmin=lb_10, xmax=ub_90), width=.2, position=position_dodge(width=.4), size=1.5) +
geom_point(pch=21, size=5, alpha=.9, position=position_dodge(width=.4)) + ga + ap + lp +
ylab("Regression coefficient \n of model parameter") +
xlab("Mean and 90% credible interval") +
scale_y_discrete(labels=rev(c(
TeX('$\\RL^{off}'),
TeX('$\\rho'),
TeX('$\\phi'),
TeX('$\\n-p_{b}'),
TeX('$\\kappa'),
TeX('$\\epsilon'),
TeX('$\\alpha_{+}'),
TeX('$\\alpha_{-}')))) +
facet_wrap(~category) + ft +
scale_fill_manual(values=c("skyblue", "blue", "orange", "red")) +
scale_x_continuous(labels=function(x) sprintf("%.1f", x)) + theme(axis.text.y = element_text(size=32)) #+ tol
all_ests_p
Spot checks of plot
# hist(learn_traces$scale.kappa.)
# hist(test_traces$scale.kappa.)
# hist(da_traces$scale.kappa.)
# hist(learn_traces$scale.rl_off.)
# hist(test_traces$scale.rl_off.)
# hist(da_traces$scale.rl_off.)
# hist(test_traces$scale.alpha_pos.)
# hist(da_traces$scale.alpha_pos.)
# hist(learn_traces$scale.epsilon.)
# hist(test_traces$scale.epsilon.)# ggsave("../paper/figs/pieces/fig8_modelpar_posteriors.png", all_ests_p ,
# width=14, height=12, dpi=500)Robustness check in univariate models the results aligned w hypotheses
summary(phi_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(phi),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.018 seconds (Warm-up)
## Chain 1: 0.034 seconds (Sampling)
## Chain 1: 0.052 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 8e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 2: 0.034 seconds (Sampling)
## Chain 2: 0.05 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.015 seconds (Warm-up)
## Chain 3: 0.032 seconds (Sampling)
## Chain 3: 0.047 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 4: 0.033 seconds (Sampling)
## Chain 4: 0.05 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(phi)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(phi) 0.0 0.1 -0.1 0.0 0.0
## sigma 1.0 0.0 0.9 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4126
## scale(phi) 0.0 1.0 4095
## sigma 0.0 1.0 3939
## mean_PPD 0.0 1.0 4023
## log-posterior 0.0 1.0 1756
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(kappa_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(kappa),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 1: 0.034 seconds (Sampling)
## Chain 1: 0.051 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 2: 0.032 seconds (Sampling)
## Chain 2: 0.049 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 8e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 3: 0.034 seconds (Sampling)
## Chain 3: 0.051 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 4: 0.033 seconds (Sampling)
## Chain 4: 0.05 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(kappa)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(kappa) -0.1 0.1 -0.1 -0.1 0.0
## sigma 1.0 0.0 0.9 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4011
## scale(kappa) 0.0 1.0 3801
## sigma 0.0 1.0 3926
## mean_PPD 0.0 1.0 3971
## log-posterior 0.0 1.0 1887
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(npb_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(not_pun_bonus),
data=m35))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 1: 0.032 seconds (Sampling)
## Chain 1: 0.049 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 2: 0.033 seconds (Sampling)
## Chain 2: 0.049 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 5e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 3: 0.033 seconds (Sampling)
## Chain 3: 0.05 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 4: 0.033 seconds (Sampling)
## Chain 4: 0.05 seconds (Total)
## Chain 4:##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(not_pun_bonus) 0.0 0.1 0.0 0.0 0.1
## sigma 1.0 0.0 0.9 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3878
## scale(not_pun_bonus) 0.0 1.0 4073
## sigma 0.0 1.0 4175
## mean_PPD 0.0 1.0 4030
## log-posterior 0.0 1.0 1984
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).Sanity check against predictions of frequentist models and check their VIF
summary(da_preds <- lm(scale(depr_anx_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## Call:
## lm(formula = scale(depr_anx_sum) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2736 -0.8104 -0.2612 0.7150 3.0819
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.124e-16 6.069e-02 0.000 1.000
## scale(kappa) -9.197e-02 7.175e-02 -1.282 0.201
## scale(alpha_pos) 5.706e-02 7.000e-02 0.815 0.416
## scale(alpha_neg) 4.365e-03 6.299e-02 0.069 0.945
## scale(phi) -7.165e-02 7.138e-02 -1.004 0.316
## scale(rho) 6.362e-03 6.506e-02 0.098 0.922
## scale(rl_off) 3.362e-03 6.608e-02 0.051 0.959
## scale(epsilon) 4.345e-02 6.302e-02 0.689 0.491
## scale(not_pun_bonus) 6.732e-02 6.285e-02 1.071 0.285
##
## Residual standard error: 1.006 on 266 degrees of freedom
## Multiple R-squared: 0.01666, Adjusted R-squared: -0.01292
## F-statistic: 0.5633 on 8 and 266 DF, p-value: 0.8076summary(rrs_preds <- lm(scale(RRS_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## Call:
## lm(formula = scale(RRS_sum) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.54304 -0.80805 -0.03143 0.61693 2.77658
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.697e-17 6.099e-02 0.000 1.000
## scale(kappa) -2.798e-03 7.211e-02 -0.039 0.969
## scale(alpha_pos) 1.444e-02 7.035e-02 0.205 0.838
## scale(alpha_neg) -2.077e-03 6.330e-02 -0.033 0.974
## scale(phi) -2.554e-02 7.173e-02 -0.356 0.722
## scale(rho) 3.520e-02 6.538e-02 0.538 0.591
## scale(rl_off) 4.736e-02 6.641e-02 0.713 0.476
## scale(epsilon) 3.551e-02 6.333e-02 0.561 0.575
## scale(not_pun_bonus) -3.289e-03 6.316e-02 -0.052 0.959
##
## Residual standard error: 1.011 on 266 degrees of freedom
## Multiple R-squared: 0.006906, Adjusted R-squared: -0.02296
## F-statistic: 0.2312 on 8 and 266 DF, p-value: 0.9849summary(learning_preds <- lm(scale(learn_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## Call:
## lm(formula = scale(learn_m) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.94486 -0.24830 0.04602 0.28117 1.60594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.684e-15 3.183e-02 0.000 1.000000
## scale(kappa) 2.549e-01 3.764e-02 6.772 8.11e-11 ***
## scale(alpha_pos) 4.301e-01 3.672e-02 11.715 < 2e-16 ***
## scale(alpha_neg) 1.557e-01 3.304e-02 4.713 3.93e-06 ***
## scale(phi) -2.687e-01 3.744e-02 -7.176 7.13e-12 ***
## scale(rho) 1.363e-01 3.412e-02 3.995 8.38e-05 ***
## scale(rl_off) -2.667e-01 3.466e-02 -7.694 2.78e-13 ***
## scale(epsilon) -3.080e-01 3.306e-02 -9.319 < 2e-16 ***
## scale(not_pun_bonus) -1.135e-01 3.297e-02 -3.444 0.000666 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5279 on 266 degrees of freedom
## Multiple R-squared: 0.7295, Adjusted R-squared: 0.7213
## F-statistic: 89.65 on 8 and 266 DF, p-value: < 2.2e-16summary(testing_preds <- lm(scale(test_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))##
## Call:
## lm(formula = scale(test_m) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0490 -0.2169 0.0521 0.3238 1.3265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.910e-15 3.284e-02 0.000 1.0000
## scale(kappa) 5.739e-02 3.883e-02 1.478 0.1406
## scale(alpha_pos) 6.386e-01 3.788e-02 16.858 < 2e-16 ***
## scale(alpha_neg) 1.718e-01 3.408e-02 5.040 8.63e-07 ***
## scale(phi) -1.834e-01 3.862e-02 -4.747 3.37e-06 ***
## scale(rho) 7.813e-02 3.520e-02 2.219 0.0273 *
## scale(rl_off) -4.875e-01 3.576e-02 -13.633 < 2e-16 ***
## scale(epsilon) -1.670e-01 3.410e-02 -4.896 1.70e-06 ***
## scale(not_pun_bonus) -3.211e-02 3.401e-02 -0.944 0.3460
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5446 on 266 degrees of freedom
## Multiple R-squared: 0.7121, Adjusted R-squared: 0.7034
## F-statistic: 82.22 on 8 and 266 DF, p-value: < 2.2e-16car::vif(learning_preds)## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519car::vif(testing_preds)## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519car::vif(da_preds)## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519car::vif(rrs_preds)## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519Comparison to choice kernel rather than RL model
m41_v1 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew34387.csv")
m41_v2 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew47450.csv")
m41 <- rbind(m41_v1, m41_v2) %>% group_by(ID) %>% slice(which.min(nll))
#write.csv(m41, "../model_res/opt/best/BEST__m41_RunHWMPRew.csv")hist(m41$learning_bias, breaks=100)
hist(m41$alpha_ck, breaks=100)
hist(m35$alpha_pos, breaks=100)
Highly correlated
ComparePars(m35$alpha_pos, m41$alpha_ck)
As expected not very correlated
ComparePars(m35$alpha_neg, m41$alpha_ck)
hist(m41$ck_off, breaks=100)
median(m41$ck_off)## [1] 0.07782407median(m35$rl_off)## [1] 0.6136962ComparePars(m35$rl_off, m41$ck_off)
Same number of pars so don’t need penalty to compare
ComparePars(m35$nll, m41$nll, "", "m35", "m41")
sum(m35$nll-m41$nll)## [1] -500.0767length(which(m35$nll < m41$nll))/length(m41$nll)## [1] 0.65818180 rather than 1/3 inits
m42_v1 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew0Inits18746.csv")
m42_v2 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew0Inits77616.csv")
m42 <- rbind(m42_v1, m42_v2) %>% group_by(ID) %>% slice(which.min(nll))Very similar but slightly worse than m41 and thus still substantially worse than m35
ComparePars(m35$nll, m42$nll, "", "m35", "m42")
sum(m35$nll-m42$nll)## [1] -510.909length(which(m35$nll < m42$nll))/length(m42$nll)## [1] 0.6618182m41_s <- read.csv("../model_res/sims/SIM_RunHWMPRew25874.csv")m41_s_learn <- m41_s %>% filter(type=="learning")
m41_s_test <- m41_s %>% filter(type=="test")pcor_ss_sim_m41 <- data.frame(m41_s_learn %>% group_by(set_size, stim_iter) %>%
summarize(m=mean(corrects), n()))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_ss_sim_m41_iter <- data.frame(m41_s_learn %>% filter(iter %in% c(1:30)) %>% group_by(stim_iter, set_size, iter) %>%
summarize(m=mean(corrects), n()))## `summarise()` has grouped output by 'stim_iter', 'set_size'. You can override
## using the `.groups` argument.sim_m41_p1 <-
ggplot(pcor_ss_sim_m41, aes(x=as.factor(stim_iter), y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
geom_jitter(data=pcor_ss_sim_m41_iter, aes(fill=as.factor(set_size)), color="black", height=0, width=.2, alpha=1, size=2, pch=21) +
geom_point(aes(fill=as.factor(set_size)), color="black", size=6, pch=21, alpha=.7) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") +
tp + ggtitle("Simulated")#ggsave("../paper/figs/pieces/supp6_sim_perf.png", sim_m41_p1, height = 3.5, width=6, dpi=300)emp_sim_perf <- emp_p1 + sim_m41_p1 emp_sim_perf
np_si_sims_summs <- m41_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_si_sims_summs_var <- m41_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter, iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))## `summarise()` has grouped output by 'stim_iter'. You can override using the
## `.groups` argument.m41_np_si_plot <- ggplot(np_si_sims_summs,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_jitter(data=np_si_sims_summs_var, aes(x=stim_iter, y=neutral_pref), height=0,
width=.2, pch=21) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, pch=21) +
annotate("text", x=3, y=.20, label="", size=8) +
annotate("text", x=7.5, y=.20, label="", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp +
ylim(-.05, .23) +
ylab("Neutral preference \n during learning") +
ggtitle("Simulated") + scale_color_gradient2() +
scale_x_continuous(breaks=seq(1, 10, 1)) + ylim(-.1, .23)## Warning in annotate("rect", xmin = 5.5, xmax = 10.5, ymin = 0, ymax = 0.22, :
## Ignoring unknown parameters: `shape`## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.m41_np_si_plot
#ggsave("../paper/figs/pieces/supp6_neut-pref_emp-and-sims.png", m41_np_si_plot, height = 4, width=11, dpi=300)Test and SI6 data
pcor_p1_test_sim_m41 <-
data.frame(m41_s_test %>% filter(phase==1) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p1_test_sim_m41_iters <-
data.frame(m41_s_test %>% filter(phase==1) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_si6_sim_m41 <-
data.frame(m41_s_learn %>% filter(stim_iter==6) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_si6_sim_m41_iters <-
data.frame(m41_s_learn %>% filter(stim_iter==6) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.pcor_p2_test_sim_m41 <-
data.frame(m41_s_test %>% filter(phase==2) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p2_test_sim_m41_iters <-
data.frame(m41_s_test %>% filter(phase==2) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.sim_p1_test_m41 <-
ggplot(pcor_p1_test_sim_m41, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p1_test_sim_m41_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(0, 1) +
tp + ggtitle("Test phase 1")
sim_si6_m41 <- ggplot(pcor_si6_sim_m41, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_si6_sim_m41_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) +
tp + ggtitle("Stimulus iteration 6")
sim_p2_test_m41 <- ggplot(pcor_p2_test_sim_m41, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p2_test_sim_m41_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) + tp +
tp + ggtitle("Test phase 2")sims_u_plot <-
sim_p1_test_m41 + sim_si6_m41 + sim_p2_test_m41 + plot_annotation(title="Simulated", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#,emps_u_plot
sims_u_plot
#ggsave("../paper/figs/pieces/supp6_simU.png", sims_u_plot, height = 4, width=11, dpi=300)hist(m35$rl_off, breaks=50)
ck_off_p <-
ggplot(m41, aes(x=ck_off)) +
geom_vline(xintercept=median(m41$ck_off), size=4) +
geom_histogram(fill="white", color="black") + ga + ap + ylab("") +
ggtitle(TeX('$\\CK^{off}')) + tp + xlab("")ck_off_p## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
median(m41$ck_off)## [1] 0.07782407median(m35$rl_off)## [1] 0.6136962Note: the learning bias par was implemented as a scalar on WM and RL such that lower values correspond to more shrinking of parameter, so plotting 1-bias so it has more intuitive interpretation as higher value = more insensitivity of WM and CK after negative PEs
lb_off_p <-
ggplot(m41, aes(x=1-learning_bias)) +
geom_vline(xintercept=median(1-m41$learning_bias), size=4) +
geom_histogram(fill="white", color="black") + ga + ap + ylab("") +
ggtitle(TeX('bias')) + tp + xlab("")
lb_off_p## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Test neutral pref
m41_test_sim_error_df <- m41_s_test %>% filter(corrects==0)
np_test_summs_m41_sim_var <- m41_test_sim_error_df %>% group_by(phase, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))## `summarise()` has grouped output by 'phase'. You can override using the
## `.groups` argument.np_test_summs_m41_sim <- m41_test_sim_error_df %>% group_by(phase) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))sim_m41_np_si_test_plot <-
ggplot(np_test_summs_m41_sim , aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_jitter(data=np_test_summs_m41_sim_var, height=0, width=.2, size=3, alpha=.8) +
geom_bar(stat="identity", color="black", alpha=.8) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ylab("") +
ggtitle("Neutral Preference at Test") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
sim_m41_np_si_test_plot
#ggsave("../paper/figs/pieces/supp6_emp_np_si.png", sim_m41_np_si_test_plot, height = 4, width=5, dpi=300)